Análisis de Sólidos Suspendidos

Autor/a
Afiliación

Ariadna Malena Seba

Fecha de publicación

8 de agosto de 2025

Resumen

Este sitio web aborda el análisis de sólidos suspendidos totales (SST) en cuerpos de agua mediante el uso de imágenes satelitales y técnicas de aprendizaje automático. Presenta una introducción teórica sobre la importancia de los SST como indicador ambiental y luego desarrolla una parte práctica con Python, donde se integran datos espectrales y mediciones reales para aplicar modelos de regresión. Se evalúa el desempeño de estos modelos y se exploran relaciones entre bandas espectrales y SST, proponiendo mejoras metodológicas para estudios futuros.

Palabras clave

GISTAQ, UTN, FRRe, Quarto, Sólidos Suspendidos

1 Sólidos suspendidos totales, por Vera Geneyer (https://github.com/VeraGeneyer)

Los sólidos suspendidos totales (TSM) son la cantidad de materia en suspensión en el agua, que incluye plancton, minerales, arena, y microorganismos. Se determinan como el residuo no filtrable de una muestra de agua. Niveles altos (TSM) pueden reducir la transparencia del agua, limitar la luz y y transportar sustancias tóxicas, afectando la vida acuática y la calidad del agua. Este parámetro, medido mediante sensores remotos, nos da información sobre el estado físico del cuerpo de agua y están relacionados con factores como la humedad, temperatura y entre otros, que es vital para detectar riesgos al ecosistema y cumplir con las normas ambientales.

1.1 Métodos tradicionales

Tabla 1: Características principales de algoritmos tradicionales para la estimación de sólidos suspendidos.
Ecuación Bandas (nm) Métricas Aguas Plataforma Referencia
\(-229.34 \left( \frac{B03}{B08} \right)^{3}+1001.65 \left( \frac{B03}{B08} \right)^{2}-1422.7 \left( \frac{B03}{B08} \right)+665.17\) B03, B08 \(R^{2}\) Embalse1 Landsat-8 [1]
\(-244.83+40.21 \cdot B01-3.67 \cdot NDWI\) B01, NDWI (B03, B08) \(R^{2}\), RMSE, d Río2 GeoEye [2]

De acuerdo a un estudio que analizó 48 cuerpos de agua, la estimación de TSM se hizo en su mayoría por modelos lineales, siendo la banda B8A la más frecuente [3].

1.2 Métodos de aprendizaje automático

El aprendizaje automático (ML) es una rama de la inteligencia artificial cuyo objetivo es desarrollar algoritmos capaces de resolver problemas mediante el análisis de datos y la creación de funciones que describen el comportamiento de fenómenos monitoreados [4]. Los modelos de aprendizaje automático más utilizados y mencionados por los investigadores para predecir la concentración de SST son:

  • Bosque Aleatorio (RF) y Refuerzo Adaptativo (AdB), modelos que se destacan por su robustez ante datos complejos y ruidosos. Estos algoritmos construyen múltiples árboles de decisión que analizan las relaciones entre características como el uso del suelo o el volumen de escorrentía y los niveles de SST [5].

  • Redes Neuronales Artificiales (ANN), copian las redes neuronales biológicas y aprenden patrones complejos en grandes volúmenes de datos, como los niveles de SST en distintas condiciones ambientales [5],

  • k-Nearest Neighbors (kNN), en sus variantes de ponderación uniforme y variable, que estima el SST en función de la cercanía en características de nuevos puntos de muestreo con datos históricos [5].

El aprendizaje automático es esencial para mejorar la precisión y rapidez en el análisis de la calidad del agua, proporcionando un monitoreo más eficiente y menos costoso en comparación con los métodos tradicionales, especialmente en áreas de difícil acceso o con datos limitados.

Tabla 2: Características principales de algoritmos de aprendizaje automático para la estimación de sólidos suspendidos.
Modelo de machine learning Software Agua Datos Métricas Referencias
Bagging y Random Forest Programa R Bahía Muestreo Prueba de normalidad multivalente Mardia-tests y Royston [4]
Regresión lineal, LASSO, regresión de vectores de soporte (SVR), K vecinos más cercanos (KNN), bosque aleatorio (RF) y redes neuronales artificiales (ANN). - Lago y embalse Sentinel-2 y UAV \(R^{2}\) [6]
Regresión lineal, regresión de vectores de soporte (SVR), K vecinos más cercanos (KNN), bosque aleatorio (RF) y redes neuronales artificiales (ANN). Programa Python Lagos Estación de monitoreo (Sensores para cada parámetro) \(R^{2}\), NSE y RMSE [5]

2 Desarrollo integral del algoritmo

En esta sección describo el proceso completo para crear un modelo que estima los sólidos suspendidos totales (SST) usando datos satelitales y de campo, desde la preparación de datos hasta la validación del algoritmo.

2.1 Procesamiento inicial de datos con Sen2Cor

Se detalla el procedimiento técnico que implementé para procesar información ambiental georreferenciada, con el objetivo de analizar el comportamiento del parámetro sólidos suspendidos (sol_sus) en una región específica del estudio (píxel 3x3).

Para esto, utilicé el lenguaje Python y la biblioteca pandas, que resulta particularmente eficiente para el manejo y transformación de datos en estructuras tabulares. Este paso fue fundamental para preparar los datos y asegurar su calidad antes del análisis.

2.1.1 Carga de datos

Primero importo la biblioteca pandas, una herramienta en Python que se utiliza para manejar datos en formato tabular (como hojas de cálculo o CSVs). Se le da el alias pd por convención, para simplificar el código.

Luego cargo dos archivos CSV con la función pd.read_csv(), la cual convierte dichos archivos en objetos del tipo DataFrame, que representan tablas en memoria, que son estructuras de datos similares a tablas (parecida a una hoja de Excel). Los conjuntos de datos cargados fueron:

  • gis_df: contiene información geográfica (latitud, longitud, pixel, etc.).
  • lab_df: contiene datos de laboratorio, incluyendo el parámetro de interés sol_sus.

Verifico la carga correcta mostrando las primeras filas con la función .head(). Es útil para ver rápidamente cómo es la estructura del archivo: qué columnas hay, qué tipo de datos, si se cargó bien.

import pandas as pd  # pandas es la biblioteca para manejar datos tabulares

# Cargar los archivos de datos
gis_df = pd.read_csv('datos/base_de_datos_gis.csv')
lab_df = pd.read_csv('datos/base_de_datos_lab.csv')

# Ver las primeras filas para asegurarse de que se cargaron bien
gis_df.head(), lab_df.head()

print("Primeras filas de gis_df:")
display(gis_df.head())

print("\nPrimeras filas de lab_df:")
display(lab_df.head())
Primeras filas de gis_df:
fecha punto pixel banda reflect longitud latitud
0 2023-05-11 1 1x1 B01 0.161744 -58.868047 -27.464687
1 2023-05-11 1 1x1 B02 0.170400 -58.868047 -27.464687
2 2023-05-11 1 1x1 B03 0.208800 -58.868047 -27.464687
3 2023-05-11 1 1x1 B04 0.251000 -58.868047 -27.464687
4 2023-05-11 1 1x1 B05 0.257656 -58.868047 -27.464687

Primeras filas de lab_df:
fecha longitud latitud param valor
0 2023-05-11 -58.868047 -27.464687 ph 6.8
1 2023-05-11 -58.868047 -27.464687 cond 141.1
2 2023-05-11 -58.868047 -27.464687 sol_sus 198.0
3 2023-05-11 -58.868047 -27.464687 turb 185.5
4 2023-05-11 -58.866729 -27.466303 ph 6.8

2.1.2 Filtrado del parámetro sólidos suspendidos (‘sol_sus’)

En el conjunto de datos del laboratorio lab_df, hay múltiples parámetros medidos (como pH, turbidez, etc.). En este caso, me interesa trabajar solamente con los datos de sólidos suspendidos, identificado como "sol_sus" en la columna param. Este filtrado selectivo lo realicé para limitar el análisis al fenómeno físico-químico de interés.

Filtré el DataFrame para quedarme solo con esas filas, y renombré la columna valor como sol_sus para que sea más claro en los siguientes pasos.

# Filtrar solo las filas donde el parámetro es 'sol_sus'
sol_sus_df = lab_df[lab_df["param"] == "sol_sus"]

# Renombrar la columna 'valor' a 'sol_sus' para que tenga sentido en el merge
sol_sus_df = sol_sus_df.rename(columns={"valor": "sol_sus"})

# Mostrar para confirmar
sol_sus_df.head()
print("Primeras filas de sol_sus_df:")
display(sol_sus_df.head())
Primeras filas de sol_sus_df:
fecha longitud latitud param sol_sus
2 2023-05-11 -58.868047 -27.464687 sol_sus 198.0
6 2023-05-11 -58.866729 -27.466303 sol_sus 150.0
10 2023-05-11 -58.864889 -27.468056 sol_sus 101.0
14 2023-05-11 -58.863367 -27.469240 sol_sus 95.0
18 2023-05-11 -58.862454 -27.470561 sol_sus 69.0

2.1.3 Reorganización de datos: pivotar bandas en columnas individuales

En este paso, convierto los valores únicos de la columna banda (como B01, B02, etc.) en nombres de columnas. Cada nueva columna contendrá los valores del parámetro reflect correspondientes a esa banda en particular. Esta operación se realiza antes de unir con los valores de sol_sus, ya que el valor de reflectancia depende de la banda, mientras que sol_sus es un dato independiente que se asignará luego por punto, fecha y ubicación.

# Pivotear la tabla para que cada banda sea una columna
gis_pivot = gis_df.pivot_table(
    index=['fecha', 'punto', 'pixel', 'latitud', 'longitud'],
    columns='banda',
    values='reflect'
).reset_index()

# Eliminar el nombre del índice de columnas generado por el pivot
gis_pivot.columns.name = None

print("Primeras filas de gis_pivot:")
display(gis_pivot.head())
Primeras filas de gis_pivot:
fecha punto pixel latitud longitud B01 B02 B03 B04 B05 B06 B07 B08 B11 B12 B8A
0 2023-05-11 1 1x1 -27.464687 -58.868047 0.161744 0.170400 0.208800 0.251000 0.257656 0.202175 0.207500 0.194000 0.118625 0.114738 0.173069
1 2023-05-11 1 3x3 -27.464687 -58.868047 0.161290 0.171256 0.206644 0.250911 0.258058 0.202170 0.207490 0.193722 0.118892 0.115462 0.173137
2 2023-05-11 2 1x1 -27.466303 -58.866729 0.149925 0.170600 0.202100 0.243200 0.247169 0.186631 0.187512 0.175800 0.114100 0.111031 0.159263
3 2023-05-11 2 3x3 -27.466303 -58.866729 0.149906 0.170378 0.202278 0.242689 0.247115 0.186739 0.188027 0.175544 0.113992 0.111210 0.159373
4 2023-05-11 3 1x1 -27.468056 -58.864889 0.149725 0.166200 0.197100 0.223800 0.218175 0.159169 0.163344 0.150700 0.114263 0.113006 0.139662

2.1.4 Unión de datos geoespaciales con datos de laboratorio

Una vez que las bandas han sido transformadas en columnas, combino esta tabla con los valores de sólidos suspendidos (sol_sus) provenientes del laboratorio. La combinación se hace usando las columnas fecha, latitud y longitud, que permiten identificar los datos correspondientes a un mismo punto geográfico y temporal.

# Realizar el merge por ubicación y fecha
df_merged = pd.merge(
    gis_pivot,
    sol_sus_df[['fecha', 'latitud', 'longitud', 'sol_sus']],
    on=['fecha', 'latitud', 'longitud'],
    how='inner'
)

print("Primeras filas del DataFrame combinado:")
display(df_merged.head())
Primeras filas del DataFrame combinado:
fecha punto pixel latitud longitud B01 B02 B03 B04 B05 B06 B07 B08 B11 B12 B8A sol_sus
0 2023-05-11 1 1x1 -27.464687 -58.868047 0.161744 0.170400 0.208800 0.251000 0.257656 0.202175 0.207500 0.194000 0.118625 0.114738 0.173069 198.0
1 2023-05-11 1 3x3 -27.464687 -58.868047 0.161290 0.171256 0.206644 0.250911 0.258058 0.202170 0.207490 0.193722 0.118892 0.115462 0.173137 198.0
2 2023-05-11 2 1x1 -27.466303 -58.866729 0.149925 0.170600 0.202100 0.243200 0.247169 0.186631 0.187512 0.175800 0.114100 0.111031 0.159263 150.0
3 2023-05-11 2 3x3 -27.466303 -58.866729 0.149906 0.170378 0.202278 0.242689 0.247115 0.186739 0.188027 0.175544 0.113992 0.111210 0.159373 150.0
4 2023-05-11 3 1x1 -27.468056 -58.864889 0.149725 0.166200 0.197100 0.223800 0.218175 0.159169 0.163344 0.150700 0.114263 0.113006 0.139662 101.0

2.1.5 Filtrado espacial por píxel de interés (3x3)

Luego de combinar los datos, aplico un filtrado adicional al DataFrame sobre la columna pixel para conservar únicamente las filas correspondientes al área geográfica designada como "3x3". Este paso reduce el dominio de análisis y permite concentrarse en una región de estudio concreta.

# Filtrar solo los datos del pixel 3x3
df_pixel_3x3 = df_merged[df_merged["pixel"] == "3x3"]

print("Primeras filas del pixel 3x3:")
display(df_pixel_3x3.head())
Primeras filas del pixel 3x3:
fecha punto pixel latitud longitud B01 B02 B03 B04 B05 B06 B07 B08 B11 B12 B8A sol_sus
1 2023-05-11 1 3x3 -27.464687 -58.868047 0.161290 0.171256 0.206644 0.250911 0.258058 0.202170 0.207490 0.193722 0.118892 0.115462 0.173137 198.0
3 2023-05-11 2 3x3 -27.466303 -58.866729 0.149906 0.170378 0.202278 0.242689 0.247115 0.186739 0.188027 0.175544 0.113992 0.111210 0.159373 150.0
5 2023-05-11 3 3x3 -27.468056 -58.864889 0.149718 0.167256 0.197322 0.223511 0.218608 0.159292 0.163228 0.151544 0.114241 0.112931 0.139783 101.0
7 2023-05-11 4 3x3 -27.469240 -58.863367 0.148441 0.163289 0.188733 0.196167 0.183854 0.139388 0.140017 0.135400 0.115602 0.113097 0.128222 95.0
9 2023-05-11 5 3x3 -27.470561 -58.862454 0.145838 0.157733 0.178767 0.173911 0.162600 0.132140 0.133229 0.126656 0.117069 0.113992 0.124548 69.0

2.1.6 Exportación de datos procesados

Finalmente, guardo el resultado como un nuevo archivo .csv dentro de la carpeta datos.

Por último, exporto el resultado a un nuevo archivo en formato .csv, mediante la función to_csv() de pandas, con el parámetro index=False para evitar que la columna de índice se incluya en el archivo de salida que pandas crea por defecto. Esto me permite utilizarlo después para visualización o análisis posterior.

# Guardar el archivo CSV dentro de la carpeta "datos"
df_pixel_3x3.to_csv('datos/datos_sol_sus_pixel_3x3.csv', index=False)

2.2 Aplicación del procesamiento a datos ACOLITE

A continuación, repetí el mismo procedimiento de procesamiento aplicado anteriormente, esta vez utilizando el archivo base_de_datos_gis_acolite.csv, generado con el procesador atmosférico ACOLITE.

ACOLITE es una herramienta diseñada específicamente para corregir efectos atmosféricos en ambientes acuáticos, a partir de imágenes satelitales. Gracias a que la estructura del archivo es similar al utilizado anteriormente (Sen2Cor), fue posible aplicar la misma lógica de filtrado, transformación y combinación para preparar los datos de reflectancia.

Código
import pandas as pd

# 1. Cargar archivos CSV 
gis_acol_df = pd.read_csv('datos/base_de_datos_gis_acolite.csv')
lab_df = pd.read_csv('datos/base_de_datos_lab.csv')

# 2. Filtrar solo el parámetro 'sol_sus' 
sol_sus_df = lab_df[lab_df["param"] == "sol_sus"]
sol_sus_df = sol_sus_df.rename(columns={"valor": "sol_sus"})

#  3. Pivotear bandas en columnas 
gis_pivot = gis_acol_df.pivot_table(
    index=['fecha', 'punto', 'latitud', 'longitud'],
    columns='banda',
    values='reflect'
).reset_index()

gis_pivot.columns.name = None  # Quitar la etiqueta 'banda'

# 4. Combinar reflectancias con sol_sus 
df_merged = pd.merge(
    gis_pivot,
    sol_sus_df[['fecha', 'latitud', 'longitud', 'sol_sus']],
    on=['fecha', 'latitud', 'longitud'],
    how='inner'
)

# 5. Guardar resultado final
df_merged.to_csv('datos/datos_sol_sus_acolite.csv', index=False)

# 6. Mostrar muestra del resultado
print("Primeras filas del DataFrame resultante:")
display(df_merged.head())
Primeras filas del DataFrame resultante:
fecha punto latitud longitud B01 B02 B03 B04 B05 B06 B07 B08 B11 B12 B8A sol_sus
0 2023-05-11 1 -27.464687 -58.868047 0.041367 0.046669 0.085630 0.128799 0.128704 0.091748 0.098365 0.083253 0.002521 0.000000 0.066612 198.0
1 2023-05-11 2 -27.466303 -58.866729 0.032235 0.047916 0.084191 0.123898 0.122481 0.079922 0.082403 0.069021 0.001669 0.000000 0.055924 150.0
2 2023-05-11 3 -27.468056 -58.864889 0.030956 0.044084 0.078444 0.104498 0.095737 0.052162 0.057475 0.044354 0.000209 0.000174 0.034767 101.0
3 2023-05-11 4 -27.469240 -58.863367 0.029646 0.040407 0.070384 0.078357 0.063937 0.032452 0.033994 0.028441 0.001279 0.000000 0.022597 95.0
4 2023-05-11 5 -27.470561 -58.862454 0.026737 0.034968 0.060701 0.056602 0.044067 0.024809 0.026684 0.019209 0.001879 0.000000 0.017930 69.0

2.3 Modelo base · Regresión lineal simple

En este análisis aplico un modelo de regresión lineal simple para estudiar la relación entre la reflectancia (banda B01) y la concentración de sólidos suspendidos (sol_sus), utilizando datos experimentales. La regresión lineal es una técnica fundamental del aprendizaje automático supervisado que permite predecir un valor continuo a partir de una o más variables independientes, y suele emplearse como modelo base frente a enfoques más complejos. A lo largo de esta sección, detallo paso a paso las acciones realizadas y explico los conceptos clave para comprender y replicar este análisis.

2.3.1 Importar librerías

En este paso, cargo las bibliotecas necesarias para procesar datos, ajustar modelos de regresión, evaluar su desempeño y visualizar los resultados.

import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt

2.3.2 Cargar datos desde un CSV

Importo el archivo .csv con los datos experimentales. Se visualizan las primeras filas para verificar que los datos se han cargado correctamente.

# Cargar el CSV
datos = pd.read_csv('datos/datos_sol_sus_acolite.csv')

# Mostrar las primeras filas para verificar
datos.head()
print("Primeras filas de datos:")
display(datos.head())
Primeras filas de datos:
fecha punto latitud longitud B01 B02 B03 B04 B05 B06 B07 B08 B11 B12 B8A sol_sus
0 2023-05-11 1 -27.464687 -58.868047 0.041367 0.046669 0.085630 0.128799 0.128704 0.091748 0.098365 0.083253 0.002521 0.000000 0.066612 198.0
1 2023-05-11 2 -27.466303 -58.866729 0.032235 0.047916 0.084191 0.123898 0.122481 0.079922 0.082403 0.069021 0.001669 0.000000 0.055924 150.0
2 2023-05-11 3 -27.468056 -58.864889 0.030956 0.044084 0.078444 0.104498 0.095737 0.052162 0.057475 0.044354 0.000209 0.000174 0.034767 101.0
3 2023-05-11 4 -27.469240 -58.863367 0.029646 0.040407 0.070384 0.078357 0.063937 0.032452 0.033994 0.028441 0.001279 0.000000 0.022597 95.0
4 2023-05-11 5 -27.470561 -58.862454 0.026737 0.034968 0.060701 0.056602 0.044067 0.024809 0.026684 0.019209 0.001879 0.000000 0.017930 69.0

2.3.3 Definir variables y particionar el conjunto

Selecciono las variables relevantes: B01 como variable independiente y sol_sus como variable dependiente. Luego divido el conjunto en dos subconjuntos: uno para entrenar el modelo y otro para probarlo, lo cual sirve para evaluar su capacidad de generalización.

# Selección de variables
X = datos[["B01"]]           # variable independiente
y = datos["sol_sus"]         # variable dependiente

# División en entrenamiento y validación
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, shuffle=False)

2.3.4 Entrenar modelo de regresión lineal

En este paso se entrena un modelo de regresión lineal usando los datos de entrenamiento. El modelo aprende la relación matemática entre la reflectancia y los sólidos suspendidos.

regressor = LinearRegression().fit(X_train, y_train)

2.3.5 Evaluar desempeño del modelo

Una vez entrenado el modelo, evaluo su desempeño usando métricas estadísticas. Estas nos permiten cuantificar qué tan bien el modelo predice los valores de sólidos suspendidos a partir de la reflectancia en los datos de prueba.

y_pred = regressor.predict(X_test)
p_rmse = mean_squared_error(y_test, y_pred)
p_r2 = r2_score(y_test, y_pred)
Métricas de desempeño
El error cuadrático medio es: 1846.79
El coeficiente de determinación (R²) es: -0.252

2.3.6 Visualizar resultados

Finalmente, se visualiza gráficamente la relación entre reflectancia y sólidos suspendidos, tanto en el conjunto de entrenamiento como en el de prueba. Esto ayuda a interpretar de forma visual cómo se ajusta el modelo a los datos reales.

fig, ax = plt.subplots(ncols=2, figsize=(10, 5), sharex=True, sharey=True)

# Gráfico entrenamiento
ax[0].plot(X_train, regressor.predict(X_train), linewidth=3, color="#17A77E", label="Modelo")
ax[0].scatter(X_train, y_train, label="Entrenamiento", color="#9D50A6", alpha=0.6)
ax[0].set(xlabel="Reflectancia", ylabel="Sol_Sus", title="Conjunto de entrenamiento")
ax[0].legend()

# Gráfico validación
ax[1].plot(X_test, y_pred, linewidth=3, color="#17A77E", label="Modelo")
ax[1].scatter(X_test, y_test, label="Validación", color="#9D50A6", alpha=0.6)
ax[1].set(xlabel="Reflectancia", ylabel="Sol_Sus", title="Conjunto de validación")
ax[1].legend()

fig.suptitle("Regresión lineal")

plt.show()

2.4 Regresión lineal individual por banda

Con el fin de profundizar el análisis, evalué la relación entre los sólidos suspendidos (sol_sus) y cada banda espectral por separado. Para esto, entrené un modelo de regresión lineal simple usando los mismos datos experimentales para cada banda.

Con este enfoque comparo la capacidad predictiva individual de cada banda mediante las métricas , ajustado y el error cuadrático medio (RMSE). Este análisis me permite identificar qué bandas tienen mejor relación lineal con los sólidos suspendidos y sentar la base para probar modelos más complejos.

Código
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
import math
from IPython.display import display

# Cargar datos
datos = pd.read_csv('datos/datos_sol_sus_acolite.csv')

# Detectar columnas de bandas
bandas = [col for col in datos.columns if col.startswith("B")]

# Lista para guardar resultados
resultados = []

# Parámetros para organización de gráficos
n_bandas = len(bandas)
ncols = 3  # Número de columnas de la grilla
nrows = math.ceil(n_bandas / ncols)  # Calculamos cuántas filas se necesitan

# Crear figura
fig, axs = plt.subplots(nrows=nrows, ncols=ncols, figsize=(12, 4 * nrows))
axs = axs.flatten()  # Asegura que podamos indexarlos como una lista

for i, banda in enumerate(bandas):
    # Variables
    X = datos[[banda]]
    y = datos["sol_sus"]

    # División entrenamiento/prueba
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, shuffle=False)

    # Ajuste del modelo
    modelo = LinearRegression().fit(X_train, y_train)
    y_train_pred = modelo.predict(X_train)

    # Métricas
    r2 = modelo.score(X_train, y_train)
    n = len(y_train)
    p = X_train.shape[1]
    r2_adj = 1 - (1 - r2) * (n - 1) / (n - p - 1)
    rmse = np.sqrt(mean_squared_error(y_train, y_train_pred))

    resultados.append({
        "Banda": banda,
        "R²": round(r2, 4),
        "R²_ajustado": round(r2_adj, 4),
        "RMSE": round(rmse, 4)
    })

    # Gráfico de entrenamiento
    ax = axs[i]
    ax.scatter(X_train, y_train, alpha=0.6, color="#9D50A6", label="Entrenamiento")
    ax.plot(X_train, y_train_pred, color="#17A77E", linewidth=1.8, label="Modelo")
    ax.set_title(f'{banda}\nR²={r2:.2f}', fontsize=10)
    ax.set_xlabel('Reflectancia', fontsize=8)
    ax.set_ylabel('Sol_Sus', fontsize=8)
    ax.tick_params(axis='both', which='major', labelsize=8)
    ax.legend(fontsize=7)
    ax.grid(True)

# Eliminar ejes sobrantes
for j in range(i + 1, len(axs)):
    fig.delaxes(axs[j])

plt.suptitle("Regresiones lineales por banda (entrenamiento)", fontsize=14)
plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()

# Tabla resumen
df_resultados = (
    pd.DataFrame(resultados)
      .round(3)
      .sort_values("R²", ascending=False)   # la mejor banda queda arriba
      .reset_index(drop=True)
)

def style_metrics(df: pd.DataFrame, caption: str):
    sty = (
        df.style
          .format(precision=3)
    )

    # Ocultar índice según versión de pandas
    if hasattr(sty, "hide"):
        sty = sty.hide(axis="index")
    else:
        sty = sty.hide_index()

    # Agregar título y estilo del título
    sty = sty.set_caption(caption)
    sty = sty.set_table_styles([{
        "selector": "caption",
        "props": [
            ("caption-side", "top"),
            ("font-weight",  "bold"),
            ("margin-bottom", "0.5em")
        ]
    }], overwrite=False)

    return sty

# Mostrar tabla  
display(style_metrics(df_resultados,
                      "Resumen de métricas por banda (entrenamiento)"))

Tabla 3: Resumen de métricas por banda (entrenamiento)
Banda R²_ajustado RMSE
B05 0.183 0.167 33.457
B06 0.180 0.164 33.523
B07 0.173 0.157 33.651
B08 0.146 0.130 34.191
B8A 0.135 0.118 34.429
B04 0.094 0.076 35.230
B12 0.038 0.019 36.303
B11 0.033 0.015 36.389
B01 0.008 -0.011 36.865
B03 0.005 -0.014 36.917
B02 0.001 -0.018 36.985
Conclusiones del análisis por banda

Los resultados muestran que ninguna banda espectral logra predecir con alta precisión los sólidos suspendidos mediante un modelo lineal simple.

Las bandas B05, B06 y B07 presentan los mejores valores de (alrededor de 0.18), lo que indica que explican menos del 20 % de la variabilidad en los datos. El resto de las bandas tienen un desempeño aún menor, con errores (RMSE) relativamente altos en comparación con los valores observados.

Esto sugiere que la relación entre reflectancia y sólidos suspendidos no es lineal en escala natural y que será necesario explorar modelos más complejos o transformaciones de datos para mejorar la predicción.

2.5 Análisis de regresión lineal con transformación logarítmica

En esta etapa aplico una transformación logarítmica natural a las variables de reflectancia y sólidos suspendidos antes de ajustar los modelos de regresión lineal. Esta transformación permite estabilizar la varianza, linealizar relaciones no lineales y reducir el impacto de valores extremos.

El procedimiento es similar al análisis anterior, pero antes de entrenar el modelo se aplica log(x) a las columnas correspondientes. Para evitar problemas con valores cero, estos se reemplazan por NaN y se excluyen del análisis. Posteriormente, entreno modelos de regresión lineal simple por banda usando las variables transformadas.

El desempeño de cada modelo se evalúa mediante métricas en escala logarítmica (R²_log, RMSE_log) y en escala original, facilitando la comparación y la interpretación de resultados.

Código
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
import math
import seaborn as sns
from IPython.display import display

# 1. ── Cargar datos ────────────────────────────────────────────────
datos = pd.read_csv("datos/datos_sol_sus_acolite.csv")
bandas = [c for c in datos.columns if c.startswith("B")]

# 2. ── Filtrar positivos (requisito para log) ──────────────────────
mask_pos   = (datos["sol_sus"] > 0) & (datos[bandas] > 0).all(axis=1)
datos_log  = datos.loc[mask_pos].copy()

# 3. ── Transformar a log-natural ──────────────────────────────────
datos_log["log_sol_sus"] = np.log(datos_log["sol_sus"])
for b in bandas:
    datos_log[f"log_{b}"] = np.log(datos_log[b])

log_bandas = [f"log_{b}" for b in bandas]

# 4. ── Train / test split (descartamos test) ───────────────────────
X_total = datos_log[log_bandas]
y_total = datos_log["log_sol_sus"]
X_train, _, y_train, _ = train_test_split(
    X_total, y_total, test_size=0.2, shuffle=False)

# 5. ── Funciones de métricas ───────────────────────────────────────
def metrics_log(model, X, y):
    y_pred = model.predict(X)
    rmse   = np.sqrt(mean_squared_error(y, y_pred))
    r2     = r2_score(y, y_pred)
    n, p   = X.shape
    r2_adj = 1 - ((1 - r2) * (n - 1)) / (n - p - 1)
    return rmse, r2, r2_adj

def metrics_original(model, X, y):
    y_pred_log  = model.predict(X).ravel()
    y_pred_orig = np.exp(y_pred_log)
    y_true_orig = np.exp(y.values)
    rmse   = np.sqrt(mean_squared_error(y_true_orig, y_pred_orig))
    r2     = r2_score(y_true_orig, y_pred_orig)
    n, p   = X.shape
    r2_adj = 1 - ((1 - r2) * (n - 1)) / (n - p - 1)
    return rmse, r2, r2_adj

# 6. ── Entrenar un modelo por banda ────────────────────────────────
resultados, resultados_nolog = [], []
lineas_log, lineas_nolog     = {}, {}

for b, lb in zip(bandas, log_bandas):
    Xtr = X_train[[lb]]
    reg = LinearRegression().fit(Xtr, y_train)

    # Métricas en log
    rmse_l, r2_l, r2a_l = metrics_log(reg, Xtr, y_train)
    resultados.append({"Banda": b, "R²_log": r2_l, 
                    "R²aj_log": r2a_l, "RMSE_log": rmse_l})

    # Métricas en original
    rmse_o, r2_o, r2a_o = metrics_original(reg, Xtr, y_train)
    resultados_nolog.append({"Banda": b,
                             "R²": r2_o, "R²aj": r2a_o, "RMSE": rmse_o})

    # Rectas ordenadas
    order        = np.argsort(Xtr.values.ravel())
    x_log_ord    = Xtr.values.ravel()[order]
    y_pred_log   = reg.predict(Xtr).ravel()[order]
    lineas_log[b] = (x_log_ord, y_pred_log)

    x_orig_ord   = np.exp(x_log_ord)
    y_pred_orig  = np.exp(y_pred_log)
    lineas_nolog[b] = (x_orig_ord, y_pred_orig)

# 7. ── DataFrames de métricas ──────────────────────────────────────
res_df_log   = (pd.DataFrame(resultados)
                  .round(3)
                  .sort_values("RMSE_log")
                  .reset_index(drop=True))

res_df_nolog = (pd.DataFrame(resultados_nolog)
                  .round(3)
                  .sort_values("RMSE")
                  .reset_index(drop=True))

# 8. ── Gráficos en escala original ────────────────────────────────
n_cols, n_rows = 3, math.ceil(len(bandas) / 3)
fig, axes = plt.subplots(n_rows, n_cols, figsize=(15, 4 * n_rows))
axes = axes.flatten()

for idx, b in enumerate(bandas):
    ax = axes[idx]
    x_orig = np.exp(X_train[f"log_{b}"])
    y_orig = np.exp(y_train)
    ax.scatter(x_orig, y_orig, alpha=0.6, color="#9D50A6", label="Entrenamiento")
    ax.plot(*lineas_nolog[b], lw=3, color="#17A77E", label="Modelo")
    ax.set(xlabel=b, ylabel="sol_sus", title=f"Banda {b} (escala original)")
    ax.legend()
    ax.grid(True)

for j in range(len(bandas), len(axes)):
    fig.delaxes(axes[j])

plt.tight_layout()
plt.show()

# 9. ── Función para estilo de tablas ───────────────────────────────
def style_metrics_sin_gradiente(df: pd.DataFrame, caption: str):
    sty = (df.style
             .format(precision=3)
             # .background_gradient(...) ← quitamos el gradiente
             .hide(axis="index")
             .set_caption(caption)
             .set_table_styles([{
                 "selector": "caption",
                 "props": [("caption-side", "top"),
                           ("font-weight",  "bold"),
                           ("margin-bottom","0.5em")]
             }], overwrite=False)
          )
    return sty


# 10 ── Mostrar tablas con el nuevo estilo pastel ───────────────────────────────
display(style_metrics_sin_gradiente(res_df_log,
        "Métricas en escala log‑log (ordenadas por RMSE_log)"))

display(style_metrics_sin_gradiente(res_df_nolog,
        "Métricas en escala original (ordenadas por RMSE)"))

Tabla 4: Métricas en escala log‑log (ordenadas por RMSE_log)
Banda R²_log R²aj_log RMSE_log
B06 0.304 0.281 0.318
B07 0.303 0.280 0.319
B05 0.272 0.248 0.326
B08 0.267 0.243 0.327
B8A 0.262 0.237 0.328
B04 0.166 0.138 0.349
B12 0.111 0.082 0.360
B11 0.065 0.034 0.369
B03 0.017 -0.016 0.378
B01 0.015 -0.018 0.379
B02 0.007 -0.026 0.380
Tabla 5: Métricas en escala original (ordenadas por RMSE)
Banda R²aj RMSE
B06 0.198 0.171 31.661
B07 0.186 0.159 31.897
B05 0.178 0.150 32.049
B8A 0.156 0.128 32.472
B08 0.151 0.123 32.567
B04 0.077 0.046 33.962
B12 0.076 0.045 33.986
B11 0.043 0.011 34.576
B01 -0.028 -0.063 35.843
B03 -0.034 -0.068 35.940
B02 -0.035 -0.070 35.969
Conclusión del ajuste log–log por banda

Los resultados muestran que al aplicar la transformación logarítmica, las bandas B06, B07 y B05 obtienen los mejores ajustes, con valores de R²_log cercanos a 0.30 y errores (RMSE_log) más bajos comparados con las demás bandas. Esto indica una mejora en la capacidad predictiva respecto al análisis en escala natural.

Sin embargo, al evaluar las métricas en escala original, el mejor desempeño alcanza un máximo cercano a 0.20, lo que aún refleja una capacidad limitada para explicar la variabilidad de los sólidos suspendidos utilizando bandas individuales.

Estos resultados sugieren que la transformación logarítmica mejora el modelo, pero que para obtener predicciones más precisas será necesario explorar combinaciones de bandas, índices espectrales o modelos más complejos.

2.6 Selección de bandas mediante regresión logarítmica y AIC

Para determinar cuáles bandas espectrales aportan más información para predecir sólidos suspendidos, ajusté modelos de regresión lineal simple en escala logarítmica para cada banda individual. Además de las métricas habituales como el coeficiente de determinación () y el error cuadrático medio (RMSE), utilicé el Criterio de Información de Akaike (AIC).

El AIC es una medida que evalúa la calidad del modelo penalizando la complejidad: un valor menor indica un mejor equilibrio entre ajuste y simplicidad. Esto ayuda a evitar sobreajuste al no favorecer modelos que mejoran el ajuste sólo por aumentar el número de parámetros sin aportar realmente información útil. Así, el AIC es una herramienta clave para seleccionar las bandas que efectivamente mejoran la predicción sin agregar complejidad innecesaria.

Código
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
import seaborn as sns
from IPython.display import display

# 1. Cargar datos y detectar bandas
datos = pd.read_csv("datos/datos_sol_sus_acolite.csv")
bandas = [c for c in datos.columns if c.startswith("B")]

# 2. Filtrar casos con valores positivos (requisito de log)
mask_pos = (datos["sol_sus"] > 0) & (datos[bandas] > 0).all(axis=1)
datos_log = datos.loc[mask_pos].copy()

# 3. Transformar a log-natural
datos_log["log_sol_sus"] = np.log(datos_log["sol_sus"])
for b in bandas:
    datos_log[f"log_{b}"] = np.log(datos_log[b])

log_bandas = [f"log_{b}" for b in bandas]

# 4. Partición entrenamiento / test 
X_total = datos_log[log_bandas]
y_total = datos_log["log_sol_sus"]

X_train, _, y_train, _ = train_test_split(
    X_total, y_total, test_size=0.20, shuffle=False
)

# 5. Funciones auxiliares (AIC + métricas)
def aic_gauss(model, X, y):
    """
    AIC para regresión lineal con ruido ~ N(0, σ²)
    AIC = n · ln(RSS/n) + 2k
    """
    n = len(y)
    k = X.shape[1] + 1  # predictores + intercepto
    rss = np.sum((y - model.predict(X)) ** 2)
    return n * np.log(rss / n) + 2 * k

def metrics_log(model, X, y):
    y_pred = model.predict(X)
    rmse = np.sqrt(mean_squared_error(y, y_pred))
    r2   = r2_score(y, y_pred)
    n, p = X.shape
    r2_adj = 1 - ((1 - r2) * (n - 1)) / (n - p - 1)
    return rmse, r2, r2_adj

# 6. Ajuste univariable banda-a-banda + AIC
resultados = []

for b in bandas:
    lb = f"log_{b}"
    Xtr = X_train[[lb]]
    reg = LinearRegression().fit(Xtr, y_train)

    rmse_log, r2_log, r2a_log = metrics_log(reg, Xtr, y_train)
    aic_val = aic_gauss(reg, Xtr, y_train)

    resultados.append({
        "Banda"    : b,
        "R²"   : r2_log,
        "R²aj" : r2a_log,
        "RMSE" : rmse_log,
        "AIC"      : aic_val
    })

# 7. Tabla ordenada por AIC
df_log  = (pd.DataFrame(resultados)
             .round(3)
             .sort_values("AIC")
             .reset_index(drop=True))

def style_metrics_sin_gradiente(df: pd.DataFrame, caption: str):
    sty = (df.style
             .format(precision=3)
             .hide(axis="index")  # oculta índice
             .set_caption(caption)
             .set_table_styles([{
                 "selector": "caption",
                 "props": [("caption-side", "top"),
                           ("font-weight",  "bold"),
                           ("margin-bottom", "0.5em")]
             }], overwrite=False)
          )
    return sty

# 9 ── Mostrar tabla con el nuevo estilo pastel AIC ───────────────────────────
display(style_metrics_sin_gradiente(
    df_log, "AIC y métricas en escala log‑log (ordenadas por AIC)"
))
Tabla 6: AIC y métricas en escala log‑log (ordenadas por AIC)
Banda R²aj RMSE AIC
B06 0.304 0.281 0.318 -69.242
B07 0.303 0.280 0.319 -69.191
B05 0.272 0.248 0.326 -67.798
B08 0.267 0.243 0.327 -67.584
B8A 0.262 0.237 0.328 -67.355
B04 0.166 0.138 0.349 -63.448
B12 0.111 0.082 0.360 -61.415
B11 0.065 0.034 0.369 -59.781
B03 0.017 -0.016 0.378 -58.184
B01 0.015 -0.018 0.379 -58.112
B02 0.007 -0.026 0.380 -57.878
Conclusión sobre la selección de bandas con AIC

Los valores obtenidos muestran que las bandas B06 y B07 son las más relevantes para predecir los sólidos suspendidos, ya que presentan los valores más bajos de AIC (-69.24 y -69.19 respectivamente), indicando el mejor equilibrio entre calidad de ajuste y simplicidad del modelo.
Estas bandas también exhiben los mejores coeficientes de determinación (R²_log alrededor de 0.30) y los menores errores (RMSE_log cerca de 0.32).

El resto de las bandas tienen valores de AIC menos favorables, lo que sugiere que, individualmente, aportan menos información útil y no justifican la complejidad del modelo. Esto refuerza la idea de que para mejorar las predicciones será necesario combinar varias bandas o explorar modelos más complejos.

2.6.1 Selección de banda espectral con validación cruzada

En esta sección selecciono la banda espectral más adecuada para estimar los sólidos suspendidos (sol_sus) usando un modelo de regresión lineal simple en escala logarítmica. A diferencia de análisis anteriores, ahora incorporo una validación cruzada que respeta la estructura temporal de los datos, para asegurar que el modelo sea confiable al aplicarse a nuevas campañas de muestreo.

¿Por qué validación Leave-One-Group-Out (LOGO)?

Los datos provienen de campañas realizadas en diferentes fechas (fecha). Esto genera grupos independientes, ya que las condiciones del río pueden variar entre campañas. Usar una validación cruzada tradicional como K-Fold o Shuffle Split mezclaría muestras de una misma campaña entre entrenamiento y prueba, lo que daría una estimación demasiado optimista del desempeño.

En cambio, LeaveOneGroupOut (LOGO) deja una fecha completa afuera en cada iteración, lo que permite evaluar si el modelo puede generalizar a campañas no vistas. Esta estrategia evita filtraciones de información (data leakage) entre los datos de entrenamiento y validación, y resulta mucho más realista en contextos ambientales, donde el objetivo es predecir nuevas campañas a partir de mediciones anteriores.

Descripción del procedimiento:

  • Preprocesamiento: aplico log() a las bandas y a la variable objetivo (sol_sus) y elimino ceros, para linearizar y estabilizar la relación.

  • Modelado: ajusto un modelo de regresión lineal simple por cada banda, normalizando los datos con StandardScaler (pipeline). Para cada modelo calculo las métricas de validación cruzada: RMSE, , R² ajustado y el AIC, que penaliza la complejidad del modelo.

  • Selección: ordeno las bandas por AIC y RMSE, y elijo la mejor.

  • Evaluación final:

    • Estimo el desempeño real sobre la campaña más reciente (test).

    • Realizo un análisis bootstrap (500 réplicas) para obtener intervalos de confianza de los parámetros y métricas.

    • Muestro gráficos de predicho vs observado en escala log y original.

Este enfoque no solo me permite elegir la banda más informativa, sino que también evalúa la estabilidad del modelo y su capacidad de generalizar en el tiempo, algo clave para el monitoreo operativo de calidad de agua.

Código
# -------------------------------------------------------------
# 0 · LIBRERÍAS
import numpy as np
import pandas as pd
from pathlib import Path
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.model_selection import LeaveOneGroupOut
from sklearn.metrics import mean_squared_error, r2_score
from IPython.display import display, Markdown

# ----------  ESTÉTICA TABLAS (sin gradiente) ------------------
def style_table(df, caption):
    return (df.style
              .format(precision=3)
              .hide(axis="index")
              .set_caption(caption)
              .set_table_styles([
                 {"selector": "caption",
                  "props": [("caption-side", "top"),
                            ("font-weight",  "bold"),
                            ("margin-bottom","0.5em"),
                            ("text-align","left")]},
                 {"selector": "th",
                  "props": [("text-align","center")]},
                 {"selector": "td",
                  "props": [("text-align","right"),
                            ("padding","0.25em 0.5em")]}
              ], overwrite=False))

# -------------------------------------------------------------
# 1 · CARGA Y PREPARACIÓN
csv_file = Path("datos/datos_sol_sus_acolite.csv")
df = pd.read_csv(csv_file).drop_duplicates()
df["fecha"] = pd.to_datetime(df["fecha"])

bandas = [c for c in df.columns if c.startswith("B")]
mask   = (df["sol_sus"]>0) & (df[bandas]>0).all(axis=1)
df     = df.loc[mask].copy()

df["log_sol_sus"] = np.log(df["sol_sus"])
for b in bandas: df[f"log_{b}"] = np.log(df[b])

log_cols = [f"log_{b}" for b in bandas]
X_total  = df[log_cols];     y_total = df["log_sol_sus"]

last_date   = df["fecha"].max()
mask_test   = df["fecha"]==last_date
mask_train  = ~mask_test
X_train,y_train = X_total[mask_train],y_total[mask_train]
X_test ,y_test  = X_total[mask_test ],y_total[mask_test ]
grupos_train    = df.loc[mask_train,"fecha"]

# -------------------------------------------------------------
# 2 · AUXILIARES
adj_r2 = lambda r2,n,p: np.nan if n-p-1<=0 else 1-(1-r2)*(n-1)/(n-p-1)
rmse   = lambda a,b: np.sqrt(mean_squared_error(a,b))
def aic_gauss(m,X,y):
    n,k=len(y),X.shape[1]+1
    rss=np.sum((y-m.predict(X))**2)
    return n*np.log(rss/n)+2*k

# -------------------------------------------------------------
# 3 · SELECCIÓN DE BANDA (AIC + LOGO)
cv = LeaveOneGroupOut()
res = []
for b in bandas:
    lb, Xb_tr = f"log_{b}", X_train[[f"log_{b}"]]
    pipe = Pipeline([("s",StandardScaler()),("m",LinearRegression())]).fit(Xb_tr,y_train)
    aic  = aic_gauss(pipe.named_steps["m"],pipe.named_steps["s"].transform(Xb_tr),y_train)
    rm,r2,ra = [],[],[]
    for tr,te in cv.split(Xb_tr,y_train,grupos_train):
        mdl = Pipeline([("s",StandardScaler()),("m",LinearRegression())]).fit(Xb_tr.iloc[tr],y_train.iloc[tr])
        pred = mdl.predict(Xb_tr.iloc[te])
        rm.append(rmse(y_train.iloc[te],pred))
        r2_ = r2_score(y_train.iloc[te],pred)
        r2.append(r2_); ra.append(adj_r2(r2_,len(te),1))
    res.append({"Banda":b,"AIC":aic,
                "RMSE_CV":np.mean(rm),"R2_CV":np.mean(r2),"AdjR2_CV":np.mean(ra)})

tbl = (pd.DataFrame(res).sort_values(["AIC","RMSE_CV"]).round(4).reset_index(drop=True))
display(style_table(tbl,"Ranking por AIC (TRAIN, log)"))

best_band = tbl.loc[0,"Banda"]
display(Markdown(f"**Mejor banda:** `{best_band}`"))

# -------------------------------------------------------------
# 4 · MODELO FINAL
lb_best   = f"log_{best_band}"
pipe_final= Pipeline([("s",StandardScaler()),("m",LinearRegression())]).fit(X_train[[lb_best]],y_train)
a,b_coef  = pipe_final.named_steps["m"].intercept_,pipe_final.named_steps["m"].coef_[0]
display(Markdown("**Ecuación log:**"));  display(Markdown(f"`log(sol_sus) = {a:.4f} + {b_coef:.4f}·log({best_band})`"))
display(Markdown("**Escala real:**"));   display(Markdown(f"`sol_sus = {np.exp(a):.4f} · {best_band}^{b_coef:.4f}`"))

# -------------------------------------------------------------
# 5 · MÉTRICAS TEST (tabla horizontal)
yhat_log  = pipe_final.predict(X_test[[lb_best]])
yhat_lin  = np.exp(yhat_log)
metrics_df = pd.DataFrame({
    "Métrica":["RMSE","R²","R² ajustado"],
    "Log":[rmse(y_test,yhat_log), r2_score(y_test,yhat_log), adj_r2(r2_score(y_test,yhat_log),len(y_test),1)],
    "Real (mg/L)":[rmse(np.exp(y_test),yhat_lin), r2_score(np.exp(y_test),yhat_lin),
                   adj_r2(r2_score(np.exp(y_test),yhat_lin),len(y_test),1)]
})
metrics_df["Log"] = metrics_df["Log"].apply(lambda x:f"{x:.3f}")
metrics_df["Real (mg/L)"]=metrics_df["Real (mg/L)"].apply(lambda x:f"{x:.3f}" if abs(x)<10 else f"{x:.2f}")
display(style_table(metrics_df,"Desempeño en TEST"))

# -------------------------------------------------------------
# 6 · BOOTSTRAP (tabla horizontal)
np.random.seed(123); n_boot=500; fechas=grupos_train.unique()
ib,cb,rm_l,r2_l,ra_l,rm_r,r2_r,ra_r=[],[],[],[],[],[],[],[]
for _ in range(n_boot):
    idx = grupos_train.isin(np.random.choice(fechas,len(fechas),replace=True))
    Xb,yb = X_train[[lb_best]].loc[idx],y_train.loc[idx]
    mdl = Pipeline([("s",StandardScaler()),("m",LinearRegression())]).fit(Xb,yb)
    ib.append(mdl.named_steps["m"].intercept_)
    cb.append(mdl.named_steps["m"].coef_[0])
    pr_log = mdl.predict(X_train[[lb_best]])
    rm_l.append(rmse(y_train,pr_log)); r2_l_ = r2_score(y_train,pr_log); r2_l.append(r2_l_)
    ra_l.append(adj_r2(r2_l_,len(y_train),1))
    pr_lin = np.exp(pr_log)
    rm_r.append(rmse(np.exp(y_train),pr_lin)); r2_r_ = r2_score(np.exp(y_train),pr_lin); r2_r.append(r2_r_)
    ra_r.append(adj_r2(r2_r_,len(y_train),1))

pct = lambda arr: np.percentile(arr,[2.5,97.5])
ci = { "inter":pct(ib), "coef":pct(cb),
       "rm_l":pct(rm_l),"r2_l":pct(r2_l),"ra_l":pct(ra_l),
       "rm_r":pct(rm_r),"r2_r":pct(r2_r),"ra_r":pct(ra_r)}

bs_df = pd.DataFrame({
    "Parámetro / Métrica":[
        "Intercepto",f"Coef {best_band}",
        "RMSE (log)","R² (log)","R² ajustado (log)",
        "RMSE (mg/L)","R² (mg/L)","R² ajustado (mg/L)"
    ],
    "Promedio":[np.mean(ib),np.mean(cb),
                np.mean(rm_l),np.mean(r2_l),np.mean(ra_l),
                np.mean(rm_r),np.mean(r2_r),np.mean(ra_r)],
    "IC 95% inf.":[ci["inter"][0],ci["coef"][0],
                   ci["rm_l"][0],ci["r2_l"][0],ci["ra_l"][0],
                   ci["rm_r"][0],ci["r2_r"][0],ci["ra_r"][0]],
    "IC 95% sup.":[ci["inter"][1],ci["coef"][1],
                   ci["rm_l"][1],ci["r2_l"][1],ci["ra_l"][1],
                   ci["rm_r"][1],ci["r2_r"][1],ci["ra_r"][1]]
})
for col in bs_df.columns[1:]:
    bs_df[col] = bs_df[col].apply(lambda x:f"{x:.4f}" if abs(float(x))<10 else f"{x:.2f}")
display(style_table(bs_df,"Bootstrap (500 réplicas) · IC 95%"))

# -------------------------------------------------------------
# 7 · GRÁFICO Predicho vs Observado  (Train + Test)
#     Izquierda → escala log   ·   Derecha → escala original
plt.rcParams.update({"axes.titlesize":10, "axes.labelsize":8,
                     "xtick.labelsize":7, "ytick.labelsize":7})

TRAIN_CLR = "#9D50A6"   # violeta
TEST_CLR  = "#FFA24B"   # naranja
ID_CLR    = "#17A77E"   # línea identidad

# --- predicciones ---
y_train_log_pred = pipe_final.predict(X_train[[lb_best]])
y_test_log_pred  = pipe_final.predict(X_test [[lb_best]])
y_train_orig, y_train_pred = np.exp(y_train), np.exp(y_train_log_pred)
y_test_orig , y_test_pred  = np.exp(y_test ), np.exp(y_test_log_pred)

fig, ax = plt.subplots(1, 2, figsize=(9, 4), sharex=False, sharey=False)

#  Escala log
ax[0].scatter(y_train, y_train_log_pred, alpha=0.6,
              color=TRAIN_CLR, edgecolor="k", label="Train (●)")
ax[0].scatter(y_test , y_test_log_pred , alpha=0.6,
              color=TEST_CLR , marker="s", edgecolor="k", label="Test  (■)")

lims = [min(y_train.min(), y_test.min()),
        max(y_train.max(), y_test.max())]
ax[0].plot(lims, lims, "--", color=ID_CLR, lw=2.5, label="y = x")
ax[0].set_xlabel("log(sol_sus) observado")
ax[0].set_ylabel("log(sol_sus) predicho")
ax[0].set_title("Escala log")
ax[0].legend(fontsize=7, frameon=False)
ax[0].grid(ls=":")

eq_log_txt = f"log(sol_sus) = {a:.3f} + {b_coef:.3f}·log({best_band})"
ax[0].text(0.04, 0.96, eq_log_txt, transform=ax[0].transAxes,
           fontsize=7, va="top",
           bbox=dict(boxstyle="round,pad=0.25", fc="white", ec="gray", lw=0.4))

#  Escala original
ax[1].scatter(y_train_orig, y_train_pred, alpha=0.6,
              color=TRAIN_CLR, edgecolor="k", label="Train (●)")
ax[1].scatter(y_test_orig , y_test_pred , alpha=0.6,
              color=TEST_CLR , marker="s", edgecolor="k", label="Test  (■)")

lims_o = [min(y_train_orig.min(), y_test_orig.min(),
              y_train_pred.min(), y_test_pred.min()),
          max(y_train_orig.max(), y_test_orig.max(),
              y_train_pred.max(), y_test_pred.max())]
ax[1].plot(lims_o, lims_o, "--", color=ID_CLR, lw=2.5, label="y = x")
ax[1].set_xlabel("sol_sus observado (mg/L)")
ax[1].set_ylabel("sol_sus predicho (mg/L)")
ax[1].set_title("Escala original")
ax[1].legend(fontsize=7, frameon=False)
ax[1].grid(ls=":")

eq_orig_txt = f"sol_sus = {np.exp(a):.4f} × {best_band}^{b_coef:.4f}"
ax[1].text(0.04, 0.96, eq_orig_txt, transform=ax[1].transAxes,
           fontsize=7, va="top",
           bbox=dict(boxstyle="round,pad=0.25", fc="white", ec="gray", lw=0.4))

fig.suptitle(f"Predicho vs Observado · Banda {best_band}", fontsize=11, y=1.05)
fig.tight_layout(pad=1.2)
plt.show()

# (opcional) restaurar tamaños por si sigues graficando después
plt.rcParams.update({"axes.titlesize":11, "axes.labelsize":9,
                     "xtick.labelsize":8, "ytick.labelsize":8})
Tabla 7: Ranking por AIC (TRAIN, log)
Banda AIC RMSE_CV R2_CV AdjR2_CV
B07 -85.218 0.291 -79.607 nan
B08 -84.907 0.285 -77.085 nan
B06 -84.048 0.298 -79.242 nan
B05 -82.491 0.325 -77.426 nan
B8A -81.287 0.310 -99.788 nan
B04 -79.479 0.340 -96.472 nan
B03 -74.454 0.341 -146.625 nan
B02 -72.724 0.344 -163.983 nan
B01 -70.995 0.342 -163.197 nan
B12 -69.509 0.332 -160.438 nan
B11 -68.415 0.317 -147.291 nan

Mejor banda: B07

Ecuación log:

log(sol_sus) = 4.2330 + 0.1938·log(B07)

Escala real:

sol_sus = 68.9254 · B07^0.1938

Tabla 8: Desempeño en TEST
Métrica Log Real (mg/L)
RMSE 0.645 70.77
-1.172 -1.326
R² ajustado -1.389 -1.558
Tabla 9: Bootstrap (500 réplicas) · IC 95%
Parámetro / Métrica Promedio IC 95% inf. IC 95% sup.
Intercepto 4.2458 4.1789 4.3614
Coef B07 0.1911 0.0600 0.2917
RMSE (log) 0.2440 0.2148 0.3839
R² (log) 0.2398 -0.7610 0.4489
R² ajustado (log) 0.2117 -0.8263 0.4285
RMSE (mg/L) 18.50 15.86 34.59
R² (mg/L) 0.1584 -1.5746 0.4585
R² ajustado (mg/L) 0.1272 -1.6700 0.4384

Conclusión del modelo con validación cruzada

Selecciono la banda B07 como la mejor candidata para estimar los sólidos suspendidos. Esta elección se basa en el valor más bajo del AIC y el menor error promedio en validación cruzada (RMSE_CV). El modelo resultante indica que, a mayor reflectancia en la banda B07, tiende a haber mayor concentración estimada de sólidos suspendidos.

Sin embargo, al evaluar el modelo con los datos más recientes (campaña de prueba), el desempeño cae considerablemente. El es negativo, lo que significa que el modelo predice peor que una estimación promedio constante. En otras palabras, el modelo no generaliza bien a nuevos datos.

Para tener una idea más completa de su estabilidad, aplico la técnica de bootstrap con 500 repeticiones. Esto muestra que si bien el coeficiente de la banda B07 es consistentemente positivo (lo que confirma una relación), el rendimiento general del modelo es bajo y muy variable.

Concluyo que, aunque B07 es la mejor banda individual, no alcanza para predecir sólidos suspendidos con precisión. Para mejorar el modelo, es necesario combinar varias bandas, usar índices espectrales o aplicar modelos más complejos.

2.7 Selección de bandas con Forward Selection basado en AIC

En esta sección realizo la selección de un conjunto de bandas espectrales para predecir la concentración de sólidos suspendidos (sol_sus) mediante un modelo de regresión lineal múltiple en escala logarítmica. Utilizo el método de Forward Selection guiado por el criterio de información de Akaike (AIC), que permite incorporar progresivamente aquellas bandas cuya inclusión mejora significativamente el ajuste del modelo sin incrementar de forma innecesaria su complejidad.

En esta sección selecciono un conjunto de bandas espectrales para predecir los sólidos suspendidos (sol_sus) mediante regresión lineal múltiple en escala logarítmica. Utilizo el método de Forward Selection, que comienza con un modelo vacío e incorpora progresivamente las bandas que más reducen el AIC, siempre que la mejora sea significativa (ΔAIC ≥ 1).

Transformo todas las variables a escala logarítmica para linealizar relaciones y reducir la variabilidad. Reservo la última fecha como conjunto de prueba, y entreno el modelo con el resto de los datos. Durante la selección, evalúo el desempeño con validación cruzada agrupada (GroupKFold), que respeta las campañas de muestreo y permite estimar la capacidad del modelo para generalizar a nuevas fechas.

Finalmente, ajusto el modelo con las bandas seleccionadas, calculo métricas en train y test, y realizo un análisis bootstrap para obtener intervalos de confianza de los coeficientes. Esto me permite construir un modelo interpretable, robusto y ajustado a las condiciones reales del monitoreo ambiental.

Código
import numpy as np
import pandas as pd
from pathlib import Path
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import GroupKFold
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt
from IPython.display import display, Markdown

# -------------------------------------------------------------
# 1 · CONFIGURACIÓN
DATOS      = Path("datos/datos_sol_sus_acolite.csv")   
N_BOOT     = 2000
DELTA_AIC  = 1.0
RNG_SEED   = 42

# -------------------------------------------------------------
# 2 · ESTILO TABLAS (robusto a versiones pandas)
def style_table(df, caption):
    sty = df.style.format(precision=3)
    if hasattr(sty, "hide"):
        sty = sty.hide(axis="index")
    elif hasattr(sty, "hide_index"):
        sty = sty.hide_index()
    # caption + estilos de alineación
    sty = (sty.set_caption(caption)
              .set_table_styles([
                  {"selector":"caption","props":[("caption-side","top"),
                                                ("font-weight","bold"),
                                                ("margin-bottom","0.5em"),
                                                ("text-align","left")]},
                  {"selector":"th","props":[("text-align","center")]},
                  {"selector":"td","props":[("text-align","right"),
                                            ("padding","0.25em 0.5em")]}
              ], overwrite=False))
    return sty

# -------------------------------------------------------------
# 3 · FUNCIONES AUXILIARES
rmse = lambda a,b: np.sqrt(mean_squared_error(a,b))
aic_g = lambda rss,n,k: n*np.log(rss/n)+2*k
def fit_metrics(Xm, ym):
    m = LinearRegression().fit(Xm, ym)
    pred = m.predict(Xm)
    n,p = Xm.shape
    rss = np.sum((ym - pred)**2)
    return {"AIC":aic_g(rss,n,p+1),
            "RMSE":rmse(ym,pred),
            "R²":r2_score(ym,pred),
            "R²aj":1-(1-r2_score(ym,pred))*(n-1)/(n-p-1),
            "mod":m}

# -------------------------------------------------------------
# 4 · CARGA Y PREPARACIÓN
df = (pd.read_csv(DATOS)
        .sort_values("fecha")
        .reset_index(drop=True))
df["fecha"] = pd.to_datetime(df["fecha"])

bandas = [c for c in df.columns if c.startswith("B")]
mask   = (df["sol_sus"]>0) & (df[bandas]>0).all(axis=1)
df     = df.loc[mask].copy()

df["log_sol_sus"] = np.log(df["sol_sus"])
for b in bandas: df[f"log_{b}"] = np.log(df[b])

X = df[[f"log_{b}" for b in bandas]]; y = df["log_sol_sus"]

last_date = df["fecha"].max()
m_test    = df["fecha"]==last_date
m_train   = ~m_test

X_tr,y_tr = X[m_train],y[m_train]
X_te,y_te = X[m_test ],y[m_test ]
groups_tr = df.loc[m_train,"fecha"]

print(f"TRAIN campañas : {groups_tr.nunique()}")
print(f"TEST  campaña  : {last_date.date()} ({m_test.sum()} muestras)")

# -------------------------------------------------------------
# 5 · FORWARD SELECTION con tabla de candidatos
remaining, selected = list(X_tr.columns), []
rss_null  = np.sum((y_tr - y_tr.mean())**2)
best_aic  = aic_g(rss_null,len(y_tr),1)
hist=[]; paso=0

while remaining:
    cand=[]
    for b in remaining:
        mets = fit_metrics(X_tr[selected+[b]].values, y_tr.values)
        cand.append((mets["AIC"], b, mets))
    cand.sort()

    df_cand = pd.DataFrame({
        "Banda":[b.replace("log_","") for _,b,_ in cand],
        "AIC":[a for a,_,_ in cand]
    })
    df_cand["ΔAIC"] = best_aic - df_cand["AIC"]
    display(Markdown(f"**Paso {paso+1} - candidatos:**"))
    display(style_table(df_cand,""))

    new_aic,best_band,best_met = cand[0]
    if best_aic - new_aic >= DELTA_AIC:
        selected.append(best_band); remaining.remove(best_band)
        best_aic=new_aic; paso+=1
        hist.append({"Paso":paso,
                     "Bandas":", ".join(c.replace("log_","") for c in selected),
                     "AIC":best_aic,
                     "RMSE":best_met["RMSE"],
                     "R²":best_met["R²"],
                     "R²aj":best_met["R²aj"]})
    else:
        break

if not selected:
    selected=[cand[0][1]]; best_aic=cand[0][0]; best_met=cand[0][2]
    hist=[{"Paso":1,"Bandas":selected[0].replace("log_",""),
           "AIC":best_aic,"RMSE":best_met["RMSE"],
           "R²":best_met["R²"],"R²aj":best_met["R²aj"]}]

hist_df=pd.DataFrame(hist)[["Paso","Bandas","AIC","RMSE","R²","R²aj"]]
display(style_table(hist_df,"Historial forward‑AIC"))

sel_log=selected; sel=[c.replace("log_","") for c in sel_log]
display(Markdown(f"**Bandas finales:** {', '.join(sel)}"))

# -------------------------------------------------------------
# 6 · GroupKFold CV
gkf=GroupKFold(n_splits=min(5,groups_tr.nunique()))
rmse_cv=[rmse(y_tr.iloc[te],
              LinearRegression().fit(X_tr.iloc[tr][sel_log],y_tr.iloc[tr])
              .predict(X_tr.iloc[te][sel_log]))
         for tr,te in gkf.split(X_tr[sel_log],y_tr,groups_tr)]
display(Markdown(f"**GroupKFold CV (log) RMSE:** {np.mean(rmse_cv):.3f} ± {np.std(rmse_cv):.3f}"))

# -------------------------------------------------------------
# 7 · MODELO FINAL + MÉTRICAS
mod = LinearRegression().fit(X_tr[sel_log], y_tr)
pred_tr_log = mod.predict(X_tr[sel_log])
pred_te_log = mod.predict(X_te[sel_log])

y_tr_lin, y_te_lin = np.exp(y_tr), np.exp(y_te)
pred_tr_lin, pred_te_lin = np.exp(pred_tr_log), np.exp(pred_te_log)

metrics_df = pd.DataFrame({
    "Métrica": ["RMSE", "R²"],
    "Train log": [rmse(y_tr, pred_tr_log), r2_score(y_tr, pred_tr_log)],
    "Test log": [rmse(y_te, pred_te_log), r2_score(y_te, pred_te_log)],
    "Train orig": [rmse(y_tr_lin, pred_tr_lin), r2_score(y_tr_lin, pred_tr_lin)],
    "Test orig": [rmse(y_te_lin, pred_te_lin), r2_score(y_te_lin, pred_te_lin)]
})
for col in metrics_df.columns[1:]:
    metrics_df[col] = metrics_df[col].apply(lambda x: f"{x:.3f}")

display(style_table(metrics_df, "Desempeño en Train y Test"))

# -------------------------------------------------------------
# 8 · BOOTSTRAP coeficientes
rng = np.random.default_rng(RNG_SEED)
dates = groups_tr.unique()
coef_bs = []

for _ in range(N_BOOT):
    idx = groups_tr.isin(rng.choice(dates, len(dates), replace=True))
    coef_bs.append(LinearRegression().fit(X_tr.loc[idx, sel_log], y_tr.loc[idx]).coef_)

coef_bs = np.array(coef_bs)
mean = coef_bs.mean(axis=0)
lo, hi = np.percentile(coef_bs, [2.5, 97.5], axis=0)

bs_df = pd.DataFrame({"Banda": sel, "Media β": mean, "IC inf": lo, "IC sup": hi})
display(style_table(bs_df, "Bootstrap (coeficientes, 95 % CI)"))

# -------------------------------------------------------------
# 9 · ECUACIONES Y BANDAS SELECCIONADAS
alpha = np.exp(mod.intercept_)
expr_log = " + ".join(f"{c:.4f}·log({b})" for c, b in zip(mod.coef_, sel))
expr_orig = " · ".join(f"{b}^{c:.4f}" for c, b in zip(mod.coef_, sel))

display(Markdown("**Bandas seleccionadas:**"))
display(Markdown(f"`{', '.join(sel)}`"))

display(Markdown("**Ecuación en escala logarítmica:**"))
display(Markdown(f"`log(sol_sus) = {mod.intercept_:.4f} + {expr_log}`"))

display(Markdown("**Ecuación en escala original:**"))
display(Markdown(f"`sol_sus = {alpha:.4f} · {expr_orig}`"))

# -------------------------------------------------------------
# 10 · GRÁFICO Predicho vs Observado (Train+Test)
plt.rcParams.update({"axes.titlesize":10,"axes.labelsize":8,
                     "xtick.labelsize":7,"ytick.labelsize":7})

TRAIN_CLR = "#9D50A6"
TEST_CLR = "#FFA24B"
ID_CLR = "#17A77E"

fig, ax = plt.subplots(1, 2, figsize=(9, 4))
# Escala log
ax[0].scatter(y_tr, pred_tr_log, alpha=0.6, color=TRAIN_CLR, edgecolor="k", label="Train (●)")
ax[0].scatter(y_te, pred_te_log, alpha=0.6, color=TEST_CLR, marker="s", edgecolor="k", label="Test  (■)")
lims = [min(y_tr.min(), y_te.min()), max(y_tr.max(), y_te.max())]
ax[0].plot(lims, lims, "--", color=ID_CLR, lw=2.5, label="y = x")
ax[0].set(xlabel="log(sol_sus) obs.", ylabel="log(sol_sus) pred.", title="Escala log")
ax[0].legend(fontsize=7, frameon=False)
ax[0].grid(ls=":")

# Escala original
ax[1].scatter(y_tr_lin, pred_tr_lin, alpha=0.6, color=TRAIN_CLR, edgecolor="k", label="Train (●)")
ax[1].scatter(y_te_lin, pred_te_lin, alpha=0.6, color=TEST_CLR, marker="s", edgecolor="k", label="Test  (■)")
lims_o = [min(y_tr_lin.min(), y_te_lin.min(), pred_tr_lin.min(), pred_te_lin.min()),
          max(y_tr_lin.max(), y_te_lin.max(), pred_tr_lin.max(), pred_te_lin.max())]
ax[1].plot(lims_o, lims_o, "--", color=ID_CLR, lw=2.5)
ax[1].set(xlabel="sol_sus obs. (mg/L)", ylabel="sol_sus pred. (mg/L)", title="Escala original")
ax[1].legend(fontsize=7, frameon=False)
ax[1].grid(ls=":")

fig.suptitle(f"Predicho vs Observado · ΔAIC ≥ {DELTA_AIC}", fontsize=11, y=1.05)
fig.tight_layout(pad=1.2)
plt.show()

# restaurar estilos si continuás graficando
plt.rcParams.update({"axes.titlesize":11,"axes.labelsize":9,
                     "xtick.labelsize":8,"ytick.labelsize":8})
TRAIN campañas : 5
TEST  campaña  : 2025-06-09 (12 muestras)

Paso 1 - candidatos:

Banda AIC ΔAIC
B07 -85.218 15.280
B08 -84.907 14.970
B06 -84.048 14.110
B05 -82.491 12.553
B8A -81.287 11.350
B04 -79.479 9.542
B03 -74.454 4.517
B02 -72.724 2.787
B01 -70.994 1.057
B12 -69.509 -0.429
B11 -68.415 -1.523

Paso 2 - candidatos:

Banda AIC ΔAIC
B02 -108.017 22.799
B01 -99.457 14.240
B03 -97.017 11.799
B11 -95.445 10.228
B12 -88.549 3.332
B8A -87.786 2.568
B06 -85.179 -0.039
B04 -84.967 -0.251
B05 -83.319 -1.899
B08 -83.234 -1.984

Paso 3 - candidatos:

Banda AIC ΔAIC
B11 -108.355 0.338
B08 -108.291 0.274
B03 -107.375 -0.641
B8A -107.364 -0.653
B04 -107.084 -0.933
B05 -106.992 -1.025
B12 -106.731 -1.286
B01 -106.713 -1.304
B06 -106.021 -1.996
Tabla 10: Historial forward‑AIC
Paso Bandas AIC RMSE R²aj
1 B07 -85.218 0.215 0.449 0.429
2 B07, B02 -108.017 0.140 0.766 0.748

Bandas finales: B07, B02

GroupKFold CV (log) RMSE: 0.125 ± 0.058

Tabla 11: Desempeño en Train y Test
Métrica Train log Test log Train orig Test orig
RMSE 0.140 0.314 9.209 36.331
0.766 0.484 0.817 0.387
Tabla 12: Bootstrap (coeficientes, 95 % CI)
Banda Media β IC inf IC sup
B07 0.635 0.429 0.713
B02 -0.698 -0.881 -0.141

Bandas seleccionadas:

B07, B02

Ecuación en escala logarítmica:

log(sol_sus) = 3.9813 + 0.6438·log(B07) + -0.7532·log(B02)

Ecuación en escala original:

sol_sus = 53.5891 · B07^0.6438 · B02^-0.7532

Conclusión del modelo seleccionado

Selecciono las bandas B07 y B02 como las mejores para estimar los sólidos suspendidos, basándome en la reducción significativa del AIC y la mejora en el error promedio de validación cruzada (RMSE_CV). El modelo indica que la reflectancia en B07 tiene un efecto positivo en la concentración estimada, mientras que B02 aporta un efecto negativo, lo que puede relacionarse con propiedades ópticas del agua y sedimentos.

Al evaluar el modelo con datos nuevos (campaña de prueba), el desempeño es razonablemente bueno, con un positivo que refleja una capacidad aceptable de generalización. El análisis bootstrap confirma la estabilidad de los coeficientes y las métricas, entregando confianza en la robustez del modelo.

En resumen, combinar estas dos bandas mejora sustancialmente la predicción en comparación con usar una sola banda, pero para lograr mayor precisión probablemente sea necesario explorar índices espectrales o modelos más avanzados.

2.8 Selección de variables mediante Forward Selection con AIC usando razones simples de bandas

En esta etapa exploro combinaciones de variables formadas por razones simples entre bandas espectrales, todas en escala logarítmica, para mejorar la predicción de los sólidos suspendidos (sol_sus). Aplico un método de Forward Selection guiado por el criterio de información de Akaike (AIC), que me permite construir un modelo parsimonioso incorporando solo aquellas razones que aportan mejoras significativas en el ajuste.

El análisis se realiza sobre el conjunto de entrenamiento definido previamente, manteniendo la partición temporal para test y respetando la estructura grupal en la validación cruzada. Además, implemento un análisis bootstrap para estimar la incertidumbre en los coeficientes y evaluar la estabilidad del modelo.

Con este procedimiento busco identificar relaciones espectrales más complejas que mejor expliquen la variabilidad de sólidos suspendidos, maximizando la capacidad predictiva sin sobreajustar, para un monitoreo ambiental más preciso y confiable.

Código
# 0 · LIBRERÍAS
import numpy as np
import pandas as pd
from pathlib import Path
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import GroupKFold
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt
from IPython.display import display, Markdown

# -------------------------------------------------------------
# 1 · CONFIGURACIÓN
DATOS     = Path("datos/datos_sol_sus_acolite.csv")     # ← ajustá si hace falta
DELTA_AIC = 2                                            # ΔAIC mínimo
N_BOOT    = 2000
RNG_SEED  = 42

# -------------------------------------------------------------
# 2 · FUNCIONES AUXILIARES
def aic_gauss(rss, n, k):           # k = nº coef + intercepto
    return n * np.log(rss / n) + 2 * k

def fit_metrics(X, y):
    """Ajusta LR y devuelve métricas + modelo."""
    m    = LinearRegression().fit(X, y)
    pred = m.predict(X)
    rss  = np.sum((y - pred) ** 2)
    n, p = X.shape
    return {
        "AIC" : aic_gauss(rss, n, p + 1),
        "RMSE": np.sqrt(mean_squared_error(y, pred)),
        "R²"  : r2_score(y, pred),
        "R²aj": 1 - (1 - r2_score(y, pred)) * (n - 1) / (n - p - 1),
        "mod" : m
    }

def rmse(a, b): return np.sqrt(mean_squared_error(a, b))

# -------------------------------------------------------------
# 3 · CARGA Y PREPARACIÓN DE DATOS
df = pd.read_csv(DATOS)
df["fecha"] = pd.to_datetime(df["fecha"])

bandas = [c for c in df.columns if c.startswith("B")]
mask   = (df["sol_sus"] > 0) & (df[bandas] > 0).all(axis=1)
df     = df.loc[mask].copy()

# variable respuesta en log
df["log_sol_sus"] = np.log(df["sol_sus"])

# -------------------------------------------------------------
# 3.1 · GENERAR RAZONES SIMPLES DE BANDAS EN LOGARITMO
ratio_cols = []
for i, b1 in enumerate(bandas):
    for b2 in bandas[i+1:]:
        col = f"log({b1}/{b2})"
        df[col] = np.log(df[b1]) - np.log(df[b2])
        ratio_cols.append(col)

X = df[ratio_cols]                # predictores = razones simples en log
y = df["log_sol_sus"]

# -------------------------------------------------------------
# 3.2 · PARTICIÓN TEMPORAL
last_date = df["fecha"].max()
m_test    = df["fecha"] == last_date
m_train   = ~m_test

X_tr, y_tr = X.loc[m_train], y.loc[m_train]
X_te, y_te = X.loc[m_test] , y.loc[m_test]
groups_tr  = df.loc[m_train, "fecha"]

print(f"TRAIN campañas : {groups_tr.nunique()}")
print(f"TEST  campaña  : {last_date.date()} ({m_test.sum()} muestras)")

# -------------------------------------------------------------
# 4 · FORWARD SELECTION (razones simples logarítmicas) con AIC + Δ
remaining = ratio_cols.copy()
selected  = []

# modelo nulo
rss_null  = np.sum((y_tr - y_tr.mean())**2)
best_aic  = aic_gauss(rss_null, len(y_tr), k=1)

hist = []
paso = 0
while remaining:
    candidatos = []
    for col in remaining:
        cols = selected + [col]
        m    = fit_metrics(X_tr[cols].values, y_tr.values)
        candidatos.append((m["AIC"], col, m))
    candidatos.sort()

    new_aic, best_col, best_met = candidatos[0]
    if best_aic - new_aic >= DELTA_AIC:     # mejora ≥ ΔAIC
        selected.append(best_col)
        remaining.remove(best_col)
        best_aic = new_aic
        paso += 1
        hist.append({
            "Paso"  : paso,
            "Bandas": ", ".join(selected),
            "AIC"   : best_aic,
            "RMSE"  : best_met["RMSE"],
            "R²"    : best_met["R²"],
            "R²aj"  : best_met["R²aj"]
        })
    else:
        break

# fallback (si nada mejora)
if not selected:
    selected = [min(candidatos, key=lambda t: t[0])[1]]
    best_aic = min(candidatos)[0]
    best_met = min(candidatos)[2]
    hist.append({
        "Paso"  : 1,
        "Bandas": selected[0],
        "AIC"   : best_aic,
        "RMSE"  : best_met["RMSE"],
        "R²"    : best_met["R²"],
        "R²aj"  : best_met["R²aj"]
    })

hist_df = (pd.DataFrame(hist)
             .loc[:, ["Paso", "Bandas", "AIC", "RMSE", "R²", "R²aj"]])

# Mostrar tabla historial con caption y sin índice
display(
    hist_df.style
    .format({"AIC": "{:.3f}", "RMSE": "{:.3f}", "R²": "{:.3f}", "R²aj": "{:.3f}"})
    .set_caption("Historial de selección")
    .hide(axis="index")
)

# Mostrar bandas seleccionadas en negrita Markdown
display(Markdown(f"**Bandas seleccionadas:**  `{', '.join(selected)}`"))

# -------------------------------------------------------------
# 5 · VALIDACIÓN GROUP‑KFOLD
gkf = GroupKFold(n_splits=min(5, groups_tr.nunique()))
rmse_cv = []
for tr_idx, te_idx in gkf.split(X_tr[selected], y_tr, groups_tr):
    m    = LinearRegression().fit(X_tr.iloc[tr_idx][selected], y_tr.iloc[tr_idx])
    pred = m.predict(X_tr.iloc[te_idx][selected])
    rmse_cv.append(rmse(y_tr.iloc[te_idx], pred))

# Mostrar GroupKFold CV en negrita Markdown
display(Markdown(f"**GroupKFold CV (log) RMSE:** {np.mean(rmse_cv):.3f} ± {np.std(rmse_cv):.3f}"))

# -------------------------------------------------------------
# 6 · MODELO FINAL Y MÉTRICAS
mod = LinearRegression().fit(X_tr[selected], y_tr)

# TRAIN
pred_tr_log = mod.predict(X_tr[selected])
rmse_tr_log = rmse(y_tr, pred_tr_log)
r2_tr_log   = r2_score(y_tr, pred_tr_log)

# TEST
pred_te_log = mod.predict(X_te[selected])
rmse_te_log = rmse(y_te, pred_te_log)
r2_te_log   = r2_score(y_te, pred_te_log)

# escala original
y_tr_lin, y_te_lin = np.exp(y_tr), np.exp(y_te)
pred_tr_lin, pred_te_lin = np.exp(pred_tr_log), np.exp(pred_te_log)

rmse_tr_lin = rmse(y_tr_lin, pred_tr_lin)
r2_tr_lin   = r2_score(y_tr_lin, pred_tr_lin)
rmse_te_lin = rmse(y_te_lin, pred_te_lin)
r2_te_lin   = r2_score(y_te_lin, pred_te_lin)

metrics_df = pd.DataFrame({
    "Métrica"    : ["RMSE", "R²"],
    "Train log"  : [rmse_tr_log, r2_tr_log],
    "Test log"   : [rmse_te_log, r2_te_log],
    "Train orig" : [rmse_tr_lin, r2_tr_lin],
    "Test orig"  : [rmse_te_lin, r2_te_lin]
})

# Mostrar tabla métricas con caption y sin índice
display(
    metrics_df.style
    .format({
        "Train log": "{:.3f}", "Test log": "{:.3f}",
        "Train orig": "{:.3f}", "Test orig": "{:.3f}"
    })
    .set_caption("Métricas finales")
    .hide(axis="index")
)

# -------------------------------------------------------------
# 7 · BOOTSTRAP AGRUPADO (coef 95 % CI)
rng    = np.random.default_rng(RNG_SEED)
dates  = groups_tr.unique()
coef_bs = []

for _ in range(N_BOOT):
    samp = rng.choice(dates, len(dates), replace=True)
    idx  = groups_tr.isin(samp)
    m_bs = LinearRegression().fit(X_tr.loc[idx, selected], y_tr.loc[idx])
    coef_bs.append(m_bs.coef_)

coef_bs = np.array(coef_bs)
mean_bs = coef_bs.mean(axis=0)
ci_bs   = np.percentile(coef_bs, [2.5, 97.5], axis=0)

bs_df = pd.DataFrame({
    "Banda"       : selected,
    "Media β"     : mean_bs,
    "IC 95% inf." : ci_bs[0],
    "IC 95% sup." : ci_bs[1]
})

# Mostrar tabla bootstrap con caption y sin índice
display(
    bs_df.style
    .format({
        "Media β": "{:.4f}",
        "IC 95% inf.": "{:.4f}",
        "IC 95% sup.": "{:.4f}"
    })
    .set_caption("Bootstrap (coeficientes, 95 % CI)")
    .hide(axis="index")
)

# -------------------------------------------------------------
# 8 · ECUACIONES
alpha      = np.exp(mod.intercept_)
expr_log   = " + ".join([f"{c:.4f}·{r}" for c, r in zip(mod.coef_, selected)])
expr_orig  = " · ".join([f"{r}^{c:.4f}" for c, r in zip(mod.coef_, selected)])

display(Markdown("Ecuación (log)"))
display(Markdown(f"`log(sol_sus) = {mod.intercept_:.4f} + {expr_log}`"))

display(Markdown("Ecuación (original)"))
display(Markdown(f"`sol_sus = {alpha:.4f} · {expr_orig}`"))

# -------------------------------------------------------------
# 9 · GRÁFICAS DE PARIDAD (TRAIN + TEST) - estilo especificado
plt.rcParams.update({"axes.titlesize":10,"axes.labelsize":8,
                     "xtick.labelsize":7,"ytick.labelsize":7})

TRAIN_CLR = "#9D50A6"
TEST_CLR = "#FFA24B"
ID_CLR = "#17A77E"

fig, ax = plt.subplots(1, 2, figsize=(9, 4))
# Escala log
ax[0].scatter(y_tr, pred_tr_log, alpha=0.6, color=TRAIN_CLR, edgecolor="k", label="Train (●)")
ax[0].scatter(y_te, pred_te_log, alpha=0.6, color=TEST_CLR, marker="s", edgecolor="k", label="Test  (■)")
lims = [min(y_tr.min(), y_te.min()), max(y_tr.max(), y_te.max())]
ax[0].plot(lims, lims, "--", color=ID_CLR, lw=2.5, label="y = x")
ax[0].set(xlabel="log(sol_sus) obs.", ylabel="log(sol_sus) pred.", title="Escala log")
ax[0].legend(fontsize=7, frameon=False)
ax[0].grid(ls=":")

# Escala original
ax[1].scatter(y_tr_lin, pred_tr_lin, alpha=0.6, color=TRAIN_CLR, edgecolor="k", label="Train (●)")
ax[1].scatter(y_te_lin, pred_te_lin, alpha=0.6, color=TEST_CLR, marker="s", edgecolor="k", label="Test  (■)")
lims_o = [min(y_tr_lin.min(), y_te_lin.min(), pred_tr_lin.min(), pred_te_lin.min()),
          max(y_tr_lin.max(), y_te_lin.max(), pred_tr_lin.max(), pred_te_lin.max())]
ax[1].plot(lims_o, lims_o, "--", color=ID_CLR, lw=2.5)
ax[1].set(xlabel="sol_sus obs. (mg/L)", ylabel="sol_sus pred. (mg/L)", title="Escala original")
ax[1].legend(fontsize=7, frameon=False)
ax[1].grid(ls=":")

fig.suptitle(f"Predicho vs Observado · ΔAIC ≥ {DELTA_AIC}", fontsize=11, y=1.05)
fig.tight_layout(pad=1.2)
plt.show()

# restaurar estilos si continuás graficando
plt.rcParams.update({"axes.titlesize":11,"axes.labelsize":9,
                     "xtick.labelsize":8,"ytick.labelsize":8})
TRAIN campañas : 5
TEST  campaña  : 2025-06-09 (12 muestras)
Tabla 13: Historial de selección
Paso Bandas AIC RMSE R²aj
1 log(B02/B07) -106.889 0.148 0.739 0.729
2 log(B02/B07), log(B05/B11) -110.321 0.135 0.784 0.767

Bandas seleccionadas: log(B02/B07), log(B05/B11)

GroupKFold CV (log) RMSE: 0.119 ± 0.048

Tabla 14: Métricas finales
Métrica Train log Test log Train orig Test orig
RMSE 0.135 0.388 9.258 45.927
0.784 0.212 0.816 0.021
Tabla 15: Bootstrap (coeficientes, 95 % CI)
Banda Media β IC 95% inf. IC 95% sup.
log(B02/B07) -0.5842 -0.6164 -0.5265
log(B05/B11) 0.0239 -0.0012 0.0348

Ecuación (log)

log(sol_sus) = 4.2142 + -0.5853·log(B02/B07) + 0.0256·log(B05/B11)

Ecuación (original)

sol_sus = 67.6397 · log(B02/B07)^-0.5853 · log(B05/B11)^0.0256

Conclusión del modelo seleccionado

En esta etapa, al usar razones simples de bandas espectrales en lugar de bandas individuales, logro un modelo con mejor ajuste en entrenamiento, con un R² de 0.784 frente a 0.766 del modelo previo. También el error promedio en la validación cruzada mejora ligeramente (RMSE CV 0.119 vs. 0.125), lo que indica una selección más eficiente de variables.

Sin embargo, al evaluar con los datos de la campaña más reciente, el desempeño baja significativamente: el R² cae a 0.212, mucho menor que el 0.484 obtenido anteriormente, y el RMSE aumenta, mostrando que la capacidad de generalización sigue siendo limitada.

El análisis bootstrap confirma que la razón log(B02/B07) tiene un coeficiente negativo estable y significativo, mientras que log(B05/B11) aporta muy poco y su coeficiente incluye cero en el intervalo de confianza, indicando que su efecto es marginal.

Desde el punto de vista físico, el modelo indica que a medida que aumenta la razón B02/B07, la concentración de sólidos suspendidos tiende a disminuir, mientras que la razón B05/B11 tiene un efecto muy débil e incierto.

A pesar de que el uso de razones simples resulta en un modelo más compacto y con una mejora modesta en el ajuste, la capacidad para predecir sólidos suspendidos en campañas futuras continúa siendo limitada.

2.8.1 Ajuste del umbral ΔAIC para selección de variables

Para evaluar la robustez y simplicidad del modelo, aumento el umbral mínimo de mejora en el criterio AIC (ΔAIC) de 2 a 4. Esto implica que la selección hacia adelante solo incorpora nuevas variables si la reducción en AIC es al menos 4, promoviendo modelos más parsimoniosos y evitando incluir variables con mejoras marginales que podrían no generalizar bien.

Este ajuste busca balancear mejor el compromiso entre ajuste y complejidad, reduciendo el riesgo de sobreajuste, especialmente relevante al aplicar el modelo a campañas de muestreo futuras.

Código
# 0 · LIBRERÍAS
import numpy as np
import pandas as pd
from pathlib import Path
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import GroupKFold
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt
from IPython.display import display, Markdown

# -------------------------------------------------------------
# 1 · CONFIGURACIÓN
DATOS     = Path("datos/datos_sol_sus_acolite.csv")     # ← ajustá si hace falta
DELTA_AIC = 4                                            # ΔAIC mínimo
N_BOOT    = 2000
RNG_SEED  = 42

# -------------------------------------------------------------
# 2 · FUNCIONES AUXILIARES
def aic_gauss(rss, n, k):           # k = nº coef + intercepto
    return n * np.log(rss / n) + 2 * k

def fit_metrics(X, y):
    """Ajusta LR y devuelve métricas + modelo."""
    m    = LinearRegression().fit(X, y)
    pred = m.predict(X)
    rss  = np.sum((y - pred) ** 2)
    n, p = X.shape
    return {
        "AIC" : aic_gauss(rss, n, p + 1),
        "RMSE": np.sqrt(mean_squared_error(y, pred)),
        "R²"  : r2_score(y, pred),
        "R²aj": 1 - (1 - r2_score(y, pred)) * (n - 1) / (n - p - 1),
        "mod" : m
    }

def rmse(a, b): return np.sqrt(mean_squared_error(a, b))

# -------------------------------------------------------------
# 3 · CARGA Y PREPARACIÓN DE DATOS
df = pd.read_csv(DATOS)
df["fecha"] = pd.to_datetime(df["fecha"])

bandas = [c for c in df.columns if c.startswith("B")]
mask   = (df["sol_sus"] > 0) & (df[bandas] > 0).all(axis=1)
df     = df.loc[mask].copy()

# variable respuesta en log
df["log_sol_sus"] = np.log(df["sol_sus"])

# -------------------------------------------------------------
# 3.1 · GENERAR RAZONES SIMPLES DE BANDAS EN LOGARITMO
ratio_cols = []
for i, b1 in enumerate(bandas):
    for b2 in bandas[i+1:]:
        col = f"log({b1}/{b2})"
        df[col] = np.log(df[b1]) - np.log(df[b2])
        ratio_cols.append(col)

X = df[ratio_cols]                # predictores = razones simples en log
y = df["log_sol_sus"]

# -------------------------------------------------------------
# 3.2 · PARTICIÓN TEMPORAL
last_date = df["fecha"].max()
m_test    = df["fecha"] == last_date
m_train   = ~m_test

X_tr, y_tr = X.loc[m_train], y.loc[m_train]
X_te, y_te = X.loc[m_test] , y.loc[m_test]
groups_tr  = df.loc[m_train, "fecha"]

print(f"TRAIN campañas : {groups_tr.nunique()}")
print(f"TEST  campaña  : {last_date.date()} ({m_test.sum()} muestras)")

# -------------------------------------------------------------
# 4 · FORWARD SELECTION (razones simples logarítmicas) con AIC + Δ
remaining = ratio_cols.copy()
selected  = []

# modelo nulo
rss_null  = np.sum((y_tr - y_tr.mean())**2)
best_aic  = aic_gauss(rss_null, len(y_tr), k=1)

hist = []
paso = 0
while remaining:
    candidatos = []
    for col in remaining:
        cols = selected + [col]
        m    = fit_metrics(X_tr[cols].values, y_tr.values)
        candidatos.append((m["AIC"], col, m))
    candidatos.sort()

    new_aic, best_col, best_met = candidatos[0]
    if best_aic - new_aic >= DELTA_AIC:     # mejora ≥ ΔAIC
        selected.append(best_col)
        remaining.remove(best_col)
        best_aic = new_aic
        paso += 1
        hist.append({
            "Paso"  : paso,
            "Bandas": ", ".join(selected),
            "AIC"   : best_aic,
            "RMSE"  : best_met["RMSE"],
            "R²"    : best_met["R²"],
            "R²aj"  : best_met["R²aj"]
        })
    else:
        break

# fallback (si nada mejora)
if not selected:
    selected = [min(candidatos, key=lambda t: t[0])[1]]
    best_aic = min(candidatos)[0]
    best_met = min(candidatos)[2]
    hist.append({
        "Paso"  : 1,
        "Bandas": selected[0],
        "AIC"   : best_aic,
        "RMSE"  : best_met["RMSE"],
        "R²"    : best_met["R²"],
        "R²aj"  : best_met["R²aj"]
    })

hist_df = (pd.DataFrame(hist)
             .loc[:, ["Paso", "Bandas", "AIC", "RMSE", "R²", "R²aj"]])

# Mostrar tabla historial con caption y sin índice
display(
    hist_df.style
    .format({"AIC": "{:.3f}", "RMSE": "{:.3f}", "R²": "{:.3f}", "R²aj": "{:.3f}"})
    .set_caption("Historial de selección")
    .hide(axis="index")
)

# Mostrar bandas seleccionadas en negrita Markdown
display(Markdown(f"**Bandas seleccionadas:**  `{', '.join(selected)}`"))

# -------------------------------------------------------------
# 5 · VALIDACIÓN GROUP‑KFOLD
gkf = GroupKFold(n_splits=min(5, groups_tr.nunique()))
rmse_cv = []
for tr_idx, te_idx in gkf.split(X_tr[selected], y_tr, groups_tr):
    m    = LinearRegression().fit(X_tr.iloc[tr_idx][selected], y_tr.iloc[tr_idx])
    pred = m.predict(X_tr.iloc[te_idx][selected])
    rmse_cv.append(rmse(y_tr.iloc[te_idx], pred))

# Mostrar GroupKFold CV en negrita Markdown
display(Markdown(f"**GroupKFold CV (log) RMSE:** {np.mean(rmse_cv):.3f} ± {np.std(rmse_cv):.3f}"))

# -------------------------------------------------------------
# 6 · MODELO FINAL Y MÉTRICAS
mod = LinearRegression().fit(X_tr[selected], y_tr)

# TRAIN
pred_tr_log = mod.predict(X_tr[selected])
rmse_tr_log = rmse(y_tr, pred_tr_log)
r2_tr_log   = r2_score(y_tr, pred_tr_log)

# TEST
pred_te_log = mod.predict(X_te[selected])
rmse_te_log = rmse(y_te, pred_te_log)
r2_te_log   = r2_score(y_te, pred_te_log)

# escala original
y_tr_lin, y_te_lin = np.exp(y_tr), np.exp(y_te)
pred_tr_lin, pred_te_lin = np.exp(pred_tr_log), np.exp(pred_te_log)

rmse_tr_lin = rmse(y_tr_lin, pred_tr_lin)
r2_tr_lin   = r2_score(y_tr_lin, pred_tr_lin)
rmse_te_lin = rmse(y_te_lin, pred_te_lin)
r2_te_lin   = r2_score(y_te_lin, pred_te_lin)

metrics_df = pd.DataFrame({
    "Métrica"    : ["RMSE", "R²"],
    "Train log"  : [rmse_tr_log, r2_tr_log],
    "Test log"   : [rmse_te_log, r2_te_log],
    "Train orig" : [rmse_tr_lin, r2_tr_lin],
    "Test orig"  : [rmse_te_lin, r2_te_lin]
})

# Mostrar tabla métricas con caption y sin índice
display(
    metrics_df.style
    .format({
        "Train log": "{:.3f}", "Test log": "{:.3f}",
        "Train orig": "{:.3f}", "Test orig": "{:.3f}"
    })
    .set_caption("Métricas finales")
    .hide(axis="index")
)

# -------------------------------------------------------------
# 7 · BOOTSTRAP AGRUPADO (coef 95 % CI)
rng    = np.random.default_rng(RNG_SEED)
dates  = groups_tr.unique()
coef_bs = []

for _ in range(N_BOOT):
    samp = rng.choice(dates, len(dates), replace=True)
    idx  = groups_tr.isin(samp)
    m_bs = LinearRegression().fit(X_tr.loc[idx, selected], y_tr.loc[idx])
    coef_bs.append(m_bs.coef_)

coef_bs = np.array(coef_bs)
mean_bs = coef_bs.mean(axis=0)
ci_bs   = np.percentile(coef_bs, [2.5, 97.5], axis=0)

bs_df = pd.DataFrame({
    "Banda"       : selected,
    "Media β"     : mean_bs,
    "IC 95% inf." : ci_bs[0],
    "IC 95% sup." : ci_bs[1]
})

# Mostrar tabla bootstrap con caption y sin índice
display(
    bs_df.style
    .format({
        "Media β": "{:.4f}",
        "IC 95% inf.": "{:.4f}",
        "IC 95% sup.": "{:.4f}"
    })
    .set_caption("Bootstrap (coeficientes, 95 % CI)")
    .hide(axis="index")
)

# -------------------------------------------------------------
# 8 · ECUACIONES
alpha      = np.exp(mod.intercept_)
expr_log   = " + ".join([f"{c:.4f}·{r}" for c, r in zip(mod.coef_, selected)])
expr_orig  = " · ".join([f"{r}^{c:.4f}" for c, r in zip(mod.coef_, selected)])

display(Markdown("Ecuación (log)"))
display(Markdown(f"`log(sol_sus) = {mod.intercept_:.4f} + {expr_log}`"))

display(Markdown("Ecuación (original)"))
display(Markdown(f"`sol_sus = {alpha:.4f} · {expr_orig}`"))

# -------------------------------------------------------------
# 9 · GRÁFICAS DE PARIDAD (TRAIN + TEST) - estilo especificado
plt.rcParams.update({"axes.titlesize":10,"axes.labelsize":8,
                     "xtick.labelsize":7,"ytick.labelsize":7})

TRAIN_CLR = "#9D50A6"
TEST_CLR = "#FFA24B"
ID_CLR = "#17A77E"

fig, ax = plt.subplots(1, 2, figsize=(9, 4))
# Escala log
ax[0].scatter(y_tr, pred_tr_log, alpha=0.6, color=TRAIN_CLR, edgecolor="k", label="Train (●)")
ax[0].scatter(y_te, pred_te_log, alpha=0.6, color=TEST_CLR, marker="s", edgecolor="k", label="Test  (■)")
lims = [min(y_tr.min(), y_te.min()), max(y_tr.max(), y_te.max())]
ax[0].plot(lims, lims, "--", color=ID_CLR, lw=2.5, label="y = x")
ax[0].set(xlabel="log(sol_sus) obs.", ylabel="log(sol_sus) pred.", title="Escala log")
ax[0].legend(fontsize=7, frameon=False)
ax[0].grid(ls=":")

# Escala original
ax[1].scatter(y_tr_lin, pred_tr_lin, alpha=0.6, color=TRAIN_CLR, edgecolor="k", label="Train (●)")
ax[1].scatter(y_te_lin, pred_te_lin, alpha=0.6, color=TEST_CLR, marker="s", edgecolor="k", label="Test  (■)")
lims_o = [min(y_tr_lin.min(), y_te_lin.min(), pred_tr_lin.min(), pred_te_lin.min()),
          max(y_tr_lin.max(), y_te_lin.max(), pred_tr_lin.max(), pred_te_lin.max())]
ax[1].plot(lims_o, lims_o, "--", color=ID_CLR, lw=2.5)
ax[1].set(xlabel="sol_sus obs. (mg/L)", ylabel="sol_sus pred. (mg/L)", title="Escala original")
ax[1].legend(fontsize=7, frameon=False)
ax[1].grid(ls=":")

fig.suptitle(f"Predicho vs Observado · ΔAIC ≥ {DELTA_AIC}", fontsize=11, y=1.05)
fig.tight_layout(pad=1.2)
plt.show()

# restaurar estilos si continuás graficando
plt.rcParams.update({"axes.titlesize":11,"axes.labelsize":9,
                     "xtick.labelsize":8,"ytick.labelsize":8})
TRAIN campañas : 5
TEST  campaña  : 2025-06-09 (12 muestras)
Tabla 16: Historial de selección
Paso Bandas AIC RMSE R²aj
1 log(B02/B07) -106.889 0.148 0.739 0.729

Bandas seleccionadas: log(B02/B07)

GroupKFold CV (log) RMSE: 0.141 ± 0.053

Tabla 17: Métricas finales
Métrica Train log Test log Train orig Test orig
RMSE 0.148 0.391 9.707 46.109
0.739 0.202 0.797 0.013
Tabla 18: Bootstrap (coeficientes, 95 % CI)
Banda Media β IC 95% inf. IC 95% sup.
log(B02/B07) -0.5622 -0.6158 -0.5001

Ecuación (log)

log(sol_sus) = 4.2877 + -0.5625·log(B02/B07)

Ecuación (original)

sol_sus = 72.8024 · log(B02/B07)^-0.5625

Impacto del aumento del umbral ΔAIC en la selección de variables

Al aumentar ΔAIC a 4, el modelo se simplifica seleccionando solo la razón log(B02/B07), manteniendo un buen ajuste en entrenamiento (R² = 0.739). Sin embargo, el desempeño en la campaña de prueba es bajo, con R² cercano a cero, indicando poca capacidad de generalización. El análisis bootstrap muestra estabilidad en el coeficiente seleccionado.

Por lo tanto, un mayor ΔAIC favorece modelos más simples, pero no mejora la predicción en nuevas campañas, por lo que es necesario explorar modelos más complejos o incluir más variables.

2.9 Selección de bandas, razones y potencias para modelar sólidos suspendidos

En esta sección amplío el conjunto de variables predictoras al incluir no solo logaritmos simples de bandas espectrales, sino también sus potencias cuadráticas y cúbicas, así como razones logarítmicas entre bandas. El objetivo es capturar relaciones no lineales y combinadas que podrían mejorar la capacidad predictiva del modelo.

Utilizo una estrategia de selección hacia adelante basada en el criterio AIC (ΔAIC ≥ 4), priorizando modelos más parsimoniosos pero informativos. La validación se realiza con partición temporal (hold-out) y GroupKFold, y se evalúa la estabilidad del modelo mediante bootstrap de coeficientes (IC del 95 %). Finalmente, se comparan las predicciones con observaciones en escalas logarítmica y original.

Código
import numpy as np
import pandas as pd
from pathlib import Path
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import GroupKFold
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt
from IPython.display import display, Markdown

# -------------------------------------------------------------
# 1 · CONFIGURACIÓN
DATOS     = Path("datos/datos_sol_sus_acolite.csv")
DELTA_AIC = 4                # umbral de mejora
N_BOOT    = 2000
RNG_SEED  = 42

# -------------------------------------------------------------
# 2 · FUNCIONES
def aic_gauss(rss, n, k):              # k = coef + intercepto
    return n*np.log(rss/n) + 2*k

def fit_metrics(X, y):
    m = LinearRegression().fit(X, y)
    pred = m.predict(X)
    rss  = np.sum((y-pred)**2)
    n,p  = X.shape
    return {
        "AIC" : aic_gauss(rss, n, p+1),
        "RMSE": np.sqrt(mean_squared_error(y, pred)),
        "R²"  : r2_score(y, pred),
        "R²aj": 1 - (1-r2_score(y,pred))*(n-1)/(n-p-1),
        "mod" : m
    }

rmse = lambda a,b: np.sqrt(mean_squared_error(a,b))

# -------------------------------------------------------------
# 3 · CARGA Y PREPARACIÓN
df = pd.read_csv(DATOS)
df["fecha"] = pd.to_datetime(df["fecha"])

bandas = [c for c in df.columns if c.startswith("B")]
df = df.loc[(df["sol_sus"]>0) & (df[bandas]>0).all(1)].copy()

df["log_sol_sus"] = np.log(df["sol_sus"])

# 3.1  log(Bx)  +  (log(Bx))²,(log(Bx))³
log_cols   = []
poly_cols  = []
for b in bandas:
    log_b = f"log_{b}"
    df[log_b] = np.log(df[b])
    df[f"{log_b}^2"] = df[log_b]**2
    df[f"{log_b}^3"] = df[log_b]**3
    log_cols.append(log_b)
    poly_cols.extend([f"{log_b}^2", f"{log_b}^3"])

# 3.2  razones en log: log(Bi/Bj) = logBi - logBj
ratio_cols = []
for i,b1 in enumerate(bandas):
    for b2 in bandas[i+1:]:
        col = f"log({b1}/{b2})"
        df[col] = df[f"log_{b1}"] - df[f"log_{b2}"]
        ratio_cols.append(col)

cand_cols = log_cols + poly_cols + ratio_cols
print(f"Total de variables candidatas: {len(cand_cols)}")

X = df[cand_cols]
y = df["log_sol_sus"]

# -------------------------------------------------------------
# 4 · HOLD‑OUT TEMPORAL
last = df["fecha"].max()
m_te  = df["fecha"]==last
m_tr  = ~m_te
X_tr, y_tr = X[m_tr], y[m_tr]
X_te, y_te = X[m_te], y[m_te]
groups_tr  = df.loc[m_tr,"fecha"]

print(f"\nTRAIN campañas : {groups_tr.nunique()}")
print(f"TEST  campaña  : {last.date()}  ({m_te.sum()} muestras)")

# -------------------------------------------------------------
# 5 · FORWARD‑AIC (Δ ≥ {DELTA_AIC})
remaining = cand_cols.copy()
selected  = []

rss0   = np.sum((y_tr - y_tr.mean())**2)
best_aic = aic_gauss(rss0, len(y_tr), 1)

hist=[]; step=0
while remaining:
    scores=[]
    for col in remaining:
        cols = selected + [col]
        mtr  = fit_metrics(X_tr[cols].values, y_tr.values)
        scores.append((mtr["AIC"], col, mtr))
    scores.sort()
    new_aic, best_col, m_info = scores[0]

    if best_aic - new_aic >= DELTA_AIC:
        selected.append(best_col)
        remaining.remove(best_col)
        best_aic = new_aic
        step += 1
        hist.append({
            "Paso":step, "Variables":", ".join(selected),
            "AIC":best_aic, "RMSE":m_info["RMSE"],
            "R²":m_info["R²"], "R²aj":m_info["R²aj"]
        })
    else:
        break

if not selected:
    best_col = min(scores, key=lambda t:t[0])[1]
    selected = [best_col]
    hist.append({"Paso":1,"Variables":best_col,
                 **{k:v for k,v in scores[0][2].items() if k!="mod"}})

hist_df = pd.DataFrame(hist)

# Mostrar tabla historial con caption y sin índice
display(
    hist_df.style
    .format({"AIC": "{:.3f}", "RMSE": "{:.3f}", "R²": "{:.3f}", "R²aj": "{:.3f}"})
    .set_caption("Historial de selección")
    .hide(axis="index")
)

# Mostrar variables seleccionadas en negrita Markdown
display(Markdown(f"**Variables finales seleccionadas:** `{', '.join(selected)}`"))

# -------------------------------------------------------------
# 6 · VALIDACIÓN GROUP‑KFOLD
gkf = GroupKFold(n_splits=min(5,groups_tr.nunique()))
rm = []
for tr_idx, te_idx in gkf.split(X_tr[selected], y_tr, groups_tr):
    m = LinearRegression().fit(X_tr.iloc[tr_idx][selected], y_tr.iloc[tr_idx])
    rm.append(rmse(y_tr.iloc[te_idx], m.predict(X_tr.iloc[te_idx][selected])))

display(Markdown(f"**GroupKFold CV (log) RMSE:** {np.mean(rm):.3f} ± {np.std(rm):.3f}"))

# -------------------------------------------------------------
# 7 · MODELO FINAL Y MÉTRICAS
mod = LinearRegression().fit(X_tr[selected], y_tr)

pred_tr = mod.predict(X_tr[selected])
pred_te = mod.predict(X_te[selected])

metrics_df = pd.DataFrame({
    "Métrica"    : ["RMSE", "R²"],
    "Train log"  : [rmse(y_tr,pred_tr), r2_score(y_tr,pred_tr)],
    "Test log"   : [rmse(y_te,pred_te), r2_score(y_te,pred_te)],
    "Train orig" : [rmse(np.exp(y_tr), np.exp(pred_tr)), r2_score(np.exp(y_tr), np.exp(pred_tr))],
    "Test orig"  : [rmse(np.exp(y_te), np.exp(pred_te)), r2_score(np.exp(y_te), np.exp(pred_te))]
})

# Mostrar tabla métricas con caption y sin índice
display(
    metrics_df.style
    .format({
        "Train log": "{:.3f}", "Test log": "{:.3f}",
        "Train orig": "{:.3f}", "Test orig": "{:.3f}"
    })
    .set_caption("Métricas finales")
    .hide(axis="index")
)

# -------------------------------------------------------------
# 8 · BOOTSTRAP AGRUPADO (IC 95 %)
rng = np.random.default_rng(RNG_SEED)
dates = groups_tr.unique()
coef_bs=[]
for _ in range(N_BOOT):
    samp = rng.choice(dates, len(dates), replace=True)
    idx  = groups_tr.isin(samp)
    m_bs = LinearRegression().fit(X_tr.loc[idx,selected], y_tr.loc[idx])
    coef_bs.append(m_bs.coef_)
coef_bs = np.asarray(coef_bs)
ci_lo, ci_hi = np.percentile(coef_bs,[2.5,97.5],axis=0)

bs_df = pd.DataFrame({
    "Variable"    : selected,
    "Media β"     : coef_bs.mean(0),
    "IC 95% inf." : ci_lo,
    "IC 95% sup." : ci_hi
})

# Mostrar tabla bootstrap con caption y sin índice
display(
    bs_df.style
    .format({
        "Media β": "{:.4f}",
        "IC 95% inf.": "{:.4f}",
        "IC 95% sup.": "{:.4f}"
    })
    .set_caption("Bootstrap (coeficientes, 95 % CI)")
    .hide(axis="index")
)

# -------------------------------------------------------------
# 9 · ECUACIÓN DEL MODELO
alpha = np.exp(mod.intercept_)
expr_log  = " + ".join([f"{c:.4f}·{v}" for c,v in zip(mod.coef_,selected)])
expr_orig = " · ".join([f"{v}^{c:.4f}" for c,v in zip(mod.coef_,selected)])

display(Markdown("Ecuación (log)"))
display(Markdown(f"`log(sol_sus) = {mod.intercept_:.4f} + {expr_log}`"))

display(Markdown("Ecuación (original)"))
display(Markdown(f"`sol_sus = {alpha:.4f} · {expr_orig}`"))

# -------------------------------------------------------------
# 10 · PARIDAD TRAIN + TEST (log y mg/L) con estilo de colores y leyenda

plt.rcParams.update({"axes.titlesize":10,"axes.labelsize":8,
                     "xtick.labelsize":7,"ytick.labelsize":7})

TRAIN_CLR = "#9D50A6"
TEST_CLR = "#FFA24B"
ID_CLR = "#17A77E"

fig, ax = plt.subplots(1, 2, figsize=(9, 4))
# Escala log
ax[0].scatter(y_tr, pred_tr, alpha=0.6, color=TRAIN_CLR, edgecolor="k", label="Train (●)")
ax[0].scatter(y_te, pred_te, alpha=0.6, color=TEST_CLR, marker="s", edgecolor="k", label="Test  (■)")
lims = [min(y_tr.min(), y_te.min()), max(y_tr.max(), y_te.max())]
ax[0].plot(lims, lims, "--", color=ID_CLR, lw=2.5, label="y = x")
ax[0].set(xlabel="log(sol_sus) obs.", ylabel="log(sol_sus) pred.", title="Escala log")
ax[0].legend(fontsize=7, frameon=False)
ax[0].grid(ls=":")

# Escala original
y_tr_lin, y_te_lin = np.exp(y_tr), np.exp(y_te)
pred_tr_lin, pred_te_lin = np.exp(pred_tr), np.exp(pred_te)

ax[1].scatter(y_tr_lin, pred_tr_lin, alpha=0.6, color=TRAIN_CLR, edgecolor="k", label="Train (●)")
ax[1].scatter(y_te_lin, pred_te_lin, alpha=0.6, color=TEST_CLR, marker="s", edgecolor="k", label="Test  (■)")
lims_o = [min(y_tr_lin.min(), y_te_lin.min(), pred_tr_lin.min(), pred_te_lin.min()),
          max(y_tr_lin.max(), y_te_lin.max(), pred_tr_lin.max(), pred_te_lin.max())]
ax[1].plot(lims_o, lims_o, "--", color=ID_CLR, lw=2.5)
ax[1].set(xlabel="sol_sus obs. (mg/L)", ylabel="sol_sus pred. (mg/L)", title="Escala original")
ax[1].legend(fontsize=7, frameon=False)
ax[1].grid(ls=":")

fig.suptitle(f"Predicho vs Observado · ΔAIC ≥ {DELTA_AIC}", fontsize=11, y=1.05)
fig.tight_layout(pad=1.2)
plt.show()

# Restaurar estilos para otras gráficas
plt.rcParams.update({"axes.titlesize":11,"axes.labelsize":9,
                     "xtick.labelsize":8,"ytick.labelsize":8})
Total de variables candidatas: 88

TRAIN campañas : 5
TEST  campaña  : 2025-06-09  (12 muestras)
Tabla 19: Historial de selección
Paso Variables AIC RMSE R²aj
1 log(B02/B07) -106.889 0.148 0.739 0.729
2 log(B02/B07), log_B8A^3 -114.331 0.126 0.812 0.797

Variables finales seleccionadas: log(B02/B07), log_B8A^3

GroupKFold CV (log) RMSE: 0.123 ± 0.034

Tabla 20: Métricas finales
Métrica Train log Test log Train orig Test orig
RMSE 0.126 0.244 8.384 33.133
0.812 0.688 0.849 0.490
Tabla 21: Bootstrap (coeficientes, 95 % CI)
Variable Media β IC 95% inf. IC 95% sup.
log(B02/B07) -0.7848 -0.9849 -0.0678
log_B8A^3 -0.0036 -0.0075 0.0085

Ecuación (log)

log(sol_sus) = 4.1224 + -0.8370·log(B02/B07) + -0.0044·log_B8A^3

Ecuación (original)

sol_sus = 61.7089 · log(B02/B07)^-0.8370 · log_B8A^3^-0.0044

2.10 Selección de variables mediante Forward Selection con AIC usando validación cruzada

Tabla 22: Selección variables (AIC, RMSE, R², R²aj)
Variables RMSE R²aj AIC
['log_ratio_B04_B05'] 0.184 0.782 0.747 -24.030
['log_ratio_B04_B05', 'log_B04'] 0.161 0.829 0.764 -24.393

Variables finales: log(B04/B05), log(B04)

Ecuación en escala logarítmica:

\(\displaystyle \log(sol\_sus) = 4.459− 0.424\,\log\left(\frac{B04}{B05}\right) − 0.116\,\log(B04)\)

Ecuación en escala original:

\(\displaystyle sol\_sus = 86.402 \times \left(\frac{B04}{B05}\right)^{-0.424} \times B04^{-0.116}\)

Tabla 23: Desempeño Train y Test
Métrica Train (log) Test (log) Train (mg/L) Test (mg/L)
RMSE 0.123 0.231 12.177 12.847
0.904 0.637 0.904 0.864
R² adj 0.897 0.547 0.897 0.829
Tabla 24: Intervalos de Confianza 95 % (Bootstrap)
Métrica Train prom Train 2.5% Train 97.5% Test prom Test 2.5% Test 97.5%
RMSE 0.115 0.083 0.144 0.171 0.062 0.244
0.908 0.847 0.954 0.677 0.170 0.960
R² adj 0.901 0.836 0.951 0.596 -0.038 0.950

2.11 Selección de variables mediante Forward Selection con AIC usando bootstrap

Tabla 25: Selección variables (Bootstrap completo)
Variables RMSE R²aj AIC
['log_ratio_B04_B05'] 0.174 0.813 0.808 -136.798
['log_ratio_B04_B05', 'log_B02'] 0.146 0.867 0.860 -148.468

Variables finales: B04/B05, B02

Ecuación en escala logarítmica:

\(\displaystyle \log(sol\_sus) = 4.459−0.371\,\log(B04/B05)−0.091\,\log(B02)\)

Ecuación en escala original:

\(\displaystyle sol\_sus = 86.402\left(\frac{B04}{B05}\right)^{-0.371} \times B02^{-0.091}\)

Tabla 26: Desempeño Train y Test
Métrica Train (log) Test (log) Train (mg/L) Test (mg/L)
RMSE 0.129 0.206 12.528 11.233
0.894 0.712 0.899 0.896
R² adj 0.886 0.640 0.891 0.870
Tabla 27: Intervalos de Confianza 95 % (Bootstrap)
Métrica Train prom Train 2.5% Train 97.5% Test prom Test 2.5% Test 97.5%
RMSE 0.120 0.089 0.147 0.161 0.065 0.221
0.897 0.821 0.945 0.713 0.191 0.950
R² adj 0.889 0.808 0.941 0.641 -0.011 0.937

2.12 Procesamiento de datos geográficos

2.12.1 Modelo gkf

\(\displaystyle \text{sol\_sus} = 61.7089 \cdot \left(\frac{B02}{B07}\right)^{-0.8370} \cdot \left(B8A^3\right)^{-0.0044}\)

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 76.45%

--- Umbral: 0 ---
Estadísticas del índice de agua:
  Min: -12.4261, Max: 17.4322, Media: 0.6814
Estadísticas de la máscara de agua:
  Valores únicos: [0 1]  (1=agua, 0=tierra, 255=nodata)
Estadísticas del parámetro sólidos suspendidos:
  Min: 0.1966, Max: 12618.6895, Media: 74.7291

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 76.45%
Umbral sugerido: 0.330

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.87%

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.74%

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.73%

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.71%

2.12.2 Modelo vc

\(\displaystyle \mathrm{sol\_sus} = 86.402 \cdot \left(\tfrac{B04}{B05}\right)^{-0.424} \cdot B04^{-0.116}\)

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 76.45%

--- Umbral: 0 ---
Estadísticas del índice de agua:
  Min: -12.4261, Max: 17.4322, Media: 0.6814
Estadísticas de la máscara de agua:
  Valores únicos: [0 1]  (1=agua, 0=tierra, 255=nodata)
Estadísticas del parámetro sólidos suspendidos:
  Min: 16.0556, Max: 2994.8086, Media: 113.4642

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.87%

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.74%

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.73%

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.71%

2.12.3 Modelo bt

\(\displaystyle \mathrm{sol\_sus} = 86.402 \cdot \left(\frac{B04}{B05}\right)^{-0.371} \cdot B02^{-0.091}\)

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 76.45%

--- Umbral: 0 ---
Estadísticas del índice de agua:
  Min: -12.4261, Max: 17.4322, Media: 0.6814
Estadísticas de la máscara de agua:
  Valores únicos: [0 1]  (1=agua, 0=tierra, 255=nodata)
Estadísticas del parámetro sólidos suspendidos:
  Min: 19.9269, Max: 51112.1016, Media: 139.5654

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.87%

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.74%

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.73%

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.71%

Referencias

[1]
D. C. R. Ramírez, «Método de Estimación de Sólidos Suspendidos Totales como Indicador de la Calidad del Agua Mediante Imágenes Satelitales». Universidad Nacional de Colombia, Facultad de Ciencias Agrarias, 2017.
[2]
G. D. J. L., J. Stephan, y Delance-Martinic., «Determinación del parámetro sólidos suspendidos totales (sst) mediante imágenes de sensores ópticos en un tramo de la cuenca media del río Bogotá (Colombia).», UD y la Geomática, pp. 19-27, 2014, doi: 10.14483/udistrital.jour.udgeo.2014.9.a02.
[3]
A. Cruz-Retana, C. Fonseca-Ortiz, R. Becerril-Piña, M. A. Gómez-Albores, M. Hernández-Téllez, y R. Arévalo-Mejia, «Characterization of spectral reflectance and TSS concentration in continental water bodies worldwide», vol. 1. pp. 4-18, 2023.
[4]
E. E. C. Vargas, «Modelamiento de Relaciones entre Parámetros Fisicoquímicos y Microbiológicos en Aguas de la Bahía Interior del Lago Titicaca-Puno (Perú) mediante Árboles de Predicción», Revista Tecnica De La Facultad De Ingenieria Universidad Del Zulia, vol. 44, pp. 154-168, ago. 2021, doi: 10.22209/rt.v44n3a02.
[5]
M. Moeini, A. Shojaeizadeh, y M. Geza, «Supervised machine learning for estimation of total suspended solids in urban watersheds», Water (Switzerland), vol. 13, ene. 2021, doi: 10.3390/w13020147.
[6]
L. S. Kupssinskü, T. T. Guimarães, E. M. D. Souza, y D. C. Zanotta, «A method for chlorophyll-a and suspended solids prediction through remote sensing and machine learning», Sensors (Switzerland), vol. 20, abr. 2020, doi: 10.3390/s20072125.

Notas

  1. Aguas lénticas.↩︎

  2. d = prueba estadística de Durbin-Watson.↩︎